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SUMMARY 

A class of similarity solutions of the hydromagnetic Navier-Stokes equations and 
Maxwell’s field equations have been found in spherical coordinates. The solutions 
may be considered as generalized Landau—Squire’s jet solutions when magneto- 
hydrodynamic interaction is present. In an attempt to soive the reduced ordinary 
differential system a perturbation expansion of the small parameter a is employed, 
where « is the ratio of kinematic viscosity to the magnetic viscosity. The physical 
meaning of the solutions is discussed and results of numerical calculation are also 
presented in tabular form. 


Symbols and notation 
Magnetic intensity. 
Electrical conductivity. 
Magnetic permeability. 
Specific heat. 
Velocity vector. 
Pressure. 
Net charge density. 
Density. 
Angle between ri, and axis of symmetry. 
Unit vectors in spherical coordinates (7, 6, ¢). 
Radial distance from origin. 
Aximuth angle. 
Magnetic viscosity = 1/(47,¢). 
Kinematic viscosity. 
Parameter defined as v/A. 
Thermal conductivity. 
Pr Prandtl number. 
I, ik-component of impulse tensor. 
f (0) Dimensionless function associated with v9. 
F(@) Dimensionless function associated with »,. 
W (0) Dimensionless function associated with Hp. 
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SUMMARY 

A class of similarity solutions of the hydromagnetic Navier-Stokes equations and 
Maxwell's field equations have been found in spherical coordinates. The solutions 
may be considered as generalized Landau-Squire’s jet solutions when magneto- 
hydrodynamic interaction is present. In an attempt to solve the reduced ordinary 
differential system a perturbation expansion of the small parameter a is employed, 
where a is the ratio of kinematic viscosity to the magnetic viscosity. The physical 
meaning of the solutions is discussed and results of numerical calculation are also 
presented in tabular form. 


Symbols and notation 
H Magnetic intensity. 
Electrical conductivity. 
Magnetic permeability. 
Specific heat. 
Velocity vector. 
Pressure. 
Net charge density. 
Density. 
Angle between ri, and axis of symmetry. 
Unit vectors in spherical coordinates (r, 6, ¢). 
Radial distance from origin. 
Aximuth angle. 
Magnetic viscosity = 1/(47,¢). 
Kinematic viscosity. 
Parameter defined as v/A. 
Thermal conductivity. 
Pr Prantl number. 
I, ik-con:ponent of impulse tensor. 
t(@) Dimensionless function associated with vp. 
F(@) Dimensionless function associated with v,. 
W (8) Dimensionless function associated with /f,. 
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2(8) Dimensionless function associated with H,. 
r'(6) Solution associated with temperature field. 
k, p Parameters determined by jet strength. 


1. Introduction 


Earzy in 1944 L. D. Landau first mentioned one exact solution of the 
Navier-Stokes equations in his book Mechanics of Compact Media [(1), 
section 19]. Physically, his solution may be interpreted as representing 
the behaviour of a submerged jet emerging from the end of a thin pipe 
into an infinite space filled with the same incompressible fluid. However, 
his work was not known to western researchers until very recently and 
from 1950 onwards Squire in England investigated this same problem 
independently and published his results in a series of papers (2, 3, 4). 

A further investigation of the same problem by V. I. Yatzeev was 
published in 1950 (5). Yatzeev first found the general form of a class of 
exact solutions of the Navier-Stokes equations and then reviewed Landau’s 
solution as a special case. Mathematically, this type of similarity solution 
is interesting because the number of exact solutions of the Navier-Stokes 
equations that can be written in closed form is extremely small. 

However, as shown in this paper, a class of exact solutions similar to 
Yatzeev’s results also exists in magnetohydrodynamics. Here, by ‘similar 
to Yatzeev’s results’, we merely mean that the velocity field may be 
expressed in the same functional form. Of course, the inclusion of hydro- 
magnetic interactions may change the physical nature of the solution 
considerably. 

The present study is mostly of mathematical interest. The physical 
implication of the solution, in general, depends upon the boundary condi- 
tions and the imposed magnetic field. However, the following discussion 
may be viewed as a generalization of the Landau—Squire solution to an 
electrically conducting fluid in the presence of an external magnetic field 
which, before hydromagnetic interaction takes place, is primarily ‘irrota- 
tional’. We postpone, for the time being, a discussion of the physical 
aspects of this problem and proceed with the analysis. 


2. Functional forms of the similarity solutions and mathematical 
reduction of the differential equations 


For an incompressible, viscous, and electrically conducting fluid, in 
steady motion, we can write the governing equations as follows: 


divv = 0, 


curl| curl vxv—e curlHx H] = vV* curl v, 
4mp 
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divH = 0, (2.3) 
] 
curl curlH—vxH]| = 0, (2.4) 
TO fle 
where H and v are magnetic field and velocity vectors and v, o, and p, 
are kinematic viscosity, conductivity, and magnetic permeability, respec- 
tively. In the derivation of (2.2), (2.3), and (2.4) we have assumed 
(i) p, & 0, ie. the net charge density is zero, 
(ii) o and pw, are constant. 


Now, adopting a right-handed spherical coordinate system (r,@,¢), we 
postulate the functional forms of the similarity solutions for v and H as 
follows: 

r r 


5) 


Bins | (Fe\(FA+ wy), (2.6) 
He]\ 7 r 


where i, and ig are unit vectors in r- and 6-directions. The functions F(@), 
f(@), Z(@), and W(@) are dimensionless. Again (2.5) and (2.6) imply that 
the motion is axisymmetric (@/@¢ = 0) and there is no azimuthal motion 
and magnetic field. 

Since the vector fields v and H are solenoidal [/2.1) and (2.3)], the 
dimensionless functions F(@), f(@), Z(@), and W(@) must satisfy the equa- 


tions : 
F(6) = —f(6)cot 0— z 


dw 
and Z(0) = — W(@)cot 6— a0 
Using (2.7) and (2.8), we may eliminate F(@) and Z(@) in the subsequent 
equations. Hence, only two unknown functions f(@) and W(@) will appear 
in the following discussion. 
Substituting (2.5), (2.6), (2.7), and (2.8) into (2.2) and integrating once, 
we obtain 


2 
(soot +55) — (w cot +) + 


dj nee 
aps sql foot 0+ i) W 5 W cot e+ a) 


_{# 
~ (de 


where c is an integration constant. 


(soot 0+ ip) +2( Foot d+ 3) +-cot 5 {foot 0+ 7a)| = 2c; (2.9) 
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Furthermore, we may rewrite (2.9) in the following form: 
wgti— if —feotd-+e]+-3.00t 05, Ago 3° if? —feotd-+e) — 
—2(55+ + 4Ww?— Af*—foot +c) =0. (2.10) 


If we consider T+ 4W*—Af*—feot0+e 


as one unknown function, then it can be obtained in closed form, as 

a cos 6—b 
sin?@ 
So far we have eles considered the equation of motion. We shall now 

take a closer look at the induction equation. Because of symmetry, the 

electrical field vanishes everywhere in the present problem (assuming no 

external field is imposed). Therefore, we have 

curl H—avxH = 4 AY cot 0+ Test 


pf +44W2—1}f?—feot0+e = (2.11) 


r? dé dé 


+a,|(foova+f)W—s(Weora+TT)li, =o, (2.12) 


or we can write 


ad +4 (Weo cot 6) = o( a s-4"), (2.13) 
where a = 4m, ov 

Consequently, we have reduced the original differential system (2.1), 
(2.2), (2.3), and (2.4) to two ordinary differential equations, (2.11) and 
(2.13). Thus, in principle, general solutions of f(@) and W(@) may be 
obtained by solving (2.11) and (2.13) simultaneously. However, solutions 
in closed form are difficult to obtain. Equation (2.11) is a generalized 
Riccati equation but its immediate solution is prevented by the presence 
of the interaction term }W*. Finally, we note that if the magnetic field 
vanishes, (2.11) reduces to the equation obtained by Landau (1) and 


Squire (2). 


3. The non-singular solutions 

It is seen that three arbitrary constants a, b, and c appear in (2.11). In 
order to limit our later discussion to those solutions which do not involve 
a singularity in the region of interest, we first investigate the relation 
between these constants and the impulse tensor. The general expression 
of the impulse tensor is given as 


“ 1 
Ty, = (p+p©+4rV. ViBu—m (+e a %— =D, E,—7 5, Hi,, 
(3.1) 
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where p® is the electromagnetic pressure. Considering a steady motion 
of incompressible and constant-property fluid with hydromagnetic inter- 
action (but without presence of electrical field), we may write the impulse 
tensor in spherical coordinates as follows: 

H? v, vgcoté 
[yg = (p+ )—20(2— : ), 


r r 


H? lv ; 
Tog = (p+ 4) + 00p—2n (; wt 7 inte 


Ht 2 9,,(2%,\ He pe 
1 (v1 pty) Be, 
lg = 9, 
Ly = 9, 


— Pe H Hp. 
4n 


Now if we postulate the pressure p(r,@) to have the functional form 
R(9) 
— 


Vg , 1 dv, 
Ip, = pvgv,—2ylr— “8 
ss il Dee 


p= vp (3.3) 
then, = expressions (2.5), (2.6), (2.7), and (2.8), we find 


ee ae (6) 4-420) 420 AT 


log = a **[R(0)+42*(0)—4WX0) +f2+2f cot 9}, 


so Bin 5 moa at +(soot e+ zn) + (5 4f + cot afi), 


dé dé 


But al Y _ prcoto+2f—2! o feote— 2s Ww (0200). (3.4) 
Clearly we have 
Ixg—Iyo = “= ( awe pf —4f*—feot é), (3.5) 
and using (2.11), we finally obtain 


a 2pv2/acosd—b — 
I44 I = | sin?@ el. 


(3.6) 


Again, if we consider the ¢-component of the momentum equation this 
takes the form 
df dR d t d 
Sa0* aot at 
Integrating once with respect to : we obtain 


aps cot 0)+Z = 0. (3.7) 


1f%(0)-+ R(0) +2 + f(@oot 0+ 4240) = d, (3.8) 
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where d is the integration constant. Multiplying (3.8) by 2pv?/r? and adding 
to (3.5), we find 


=| RO)+ 42 8) +48) + 27 | ss | Sag +4]: (3.9) 
Thus 


lag = S| R (8) + +9240) +4W 40/4255] — Bel eee a —e+4], 


r? sin?0 
(3.10) 
Similarly ‘rom (3.4), (3.8), and (3.9), we may obtain 
Ing = © [R(O)—4W%(0)+4.24(0)+f?-+ 2f cot 6] 

SEER acosd—b d 

~ gin? erie: 
CTex 
Ox}, 


Now, from the condition zn 6 


we obtain the relation 


cot 6 l 


olay — Ig + rand ee ln 


) 
= aang arsin Blog = 0. (3.12) 


Since we know from (3.4) that r°Jp, is a function of @ only, 


1 


2sin 0I,.) = 0 
aap a!" mn Oly) 


and hence 


i as acos@—b Pipl se acosO—b — s 
I, = a |(“ snag —¢+d)oot + — 76 do? re oe d 


(3.13) 


r? sin 6 


-(* + 2¢ cos “) 


Now, on referring back to (3.8), we notice that the integration constant 
d must depend on the constants a, b, and c. This is due to the fact that 
the differential equation governing the velocity field f(@) is of third order. 
In order to show that a, b, c, and d are not independent, we will use the 
condition al, 


—H == Q, (3.14) 
OX, 


Thus 


1 ‘ 1 s 
mer +> = 0, Le 
rsin6 <¢r sin 6L,,) + sin 6 S(rsin ly) ++ ~ (loo 144) = (3.15) 


or 


3 oe ; 
rsind é p(rsin lg) = = —— — — (a+20086) = +5 (2c). 


Teo Tyg = r2 sin @ d@ 
(3.16) 





MAGNETOHYDRODYNAMIC NAVIER-STOKES EQUATIONS 
But from (3.10) and (3.11), we have 
y2 
lotIyg = > (2d). 


Comparing; (3.16), (3.17), we conclude that 
c= d. 
The results (3.10), (3.11), and (3.12) are very useful in the general 
discussion of the present study. It can be seen immediately that only a 
limited number of selections of the arbitrary constants a, b, and ¢ will 
provide solutions of physical interest. For simplicity, a brief discussion 
is presented in tabular form as follows: 


(3.17) 


(3.18) 





Case (1) 
a=b=c=0 


Case (2) 
a=b = —2% 


Case (3) 
a=c=0,b0 





0 


Regular for 


—n7<@O0< a2 


Ig = 


for §@= 0 


2pv*c 


r2 


Regular for 
—7r< 8< 0 
0O<O0< 97 
Singular at 6 = 0 
and @ = =F 





Regular for 


—@F < 


I, = 
at 6 = 


Regular for 
—7<6<0 
0O<0<-f 
Singular at 6 = 0 
and @ = tT 





Regular for 


ne 
Ig sar 
at@=—0 


0 





Remarks 





Landau and Squire's 
(1) (2) submerged 
jet for non-con- 
ducting fluid be- 
longs to this class 





Squire's (3) wall jet 
solution belongs to 
this class 





Squire’s (4) radial jet 
solution belongs to 
this class in which 
along @=0, the 
singularity appears 





In this paper the discussion will be restricted to Cases (1) and (2) for 
which no singularity occurs in the region of interest. 


4. Method of solution 

In section 2 the possibility of existence of similarity solutions postulated 
in (2.5) and (2.6) has been verified. As a consequence, the governing 
differential equations reduce to the form 


Gf _ ses ims acos@—b 
6 $f?--3W +f cot 0+ ————¢, 
@G@w id 


oe ee 
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It is seen that general solutions in closed form are not obtainable because 
of the coupling through the non-linear term }W? in the first equation. 
For this reason, from now on we will restrict our discussions to the special 
situation «a < 1. This dimensionless parameter is the ratio of the kinematic 
viscosity to the magnetic viscosity. In fact, a is fairly small in many 
aerophysical problems, its order of magnitude being about 10-5 to 10-+. 
Therefore, the situation « < 1 is quite realistic unless the fluid is highly 
conducting and viscous. 


To study the case « < 1, we shall use the following perturbation expan- 
eR $(8) = fol) +0f(0)+...+0"f (8), 
W(0) = W(0)+aW,(0)+...+0"W, (8). (4.1) 
As a result of direct substitution of (4.1) into (2.11) and (2.13), a family 
of subsidiary Pe. are obtained, thus 
d? Mo 
“ae 


. a 4 dW, wo 


AM cot #) = 0, (4.2) 


lle (W, cot #) = asad HN aM _wto_y 7 (4.4) 


na fo dé ‘hae Td 


df, acos6—b c+afi+ _4pe . 
Bae OTM ee wie 


df, 
d@ 


140 Wo a9’ 


= (f,+cot 0) f,—W, W,, 


Us — (f+ cot 6) fat 3/1—(4W3+M M). (4.7) 


Before proceeding to solve the above set of equations, we shall first 
examine the zero-order solutions W,(@) and f,(@) so that the reasonable- 
ness of this expansion procedure may be verified. 

It is seen that W,(@) may be obtained readily from (4.2) by direct 
integrations. Thus A cos 6— B 


~ gin@ 


w= (4.8) 


where A, B are arbitrary constants. 

The physical significance of W,(@), according to (4.1), is associated with 
the @-component of magnetic field in case « = 0. In the present problem, 
a = 0 implies that the fluid is non-conducting since the inviscid case has 
been excluded automatically by the introduction of the expressions (2.5) 
and (2.6). With this in mind, one would expect that there should be no 
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interaction between the flow field and the magnetic field. Hence, in (4.5), 
the term 4W? should be absorbable into the term 


acos 6—b 
sin?@ ” 
since a and 6 are arbitrary constants. 

In fact, this can be verified immediately from (4.8). Therefore, there 
is no loss of generality in dropping the term }W% in (4.5) so long as the 
constants a, b, and c are left arbitrary. 

The corresponding magnetic field may be written as 


a. | Se bee 
cghtis wer] r sin? He} 


name(s) "4 2) “ 
Me? rN \ Me 


which is irrotational and produces no ponderomotive force. In the follow- 
ing discussion, we require the behaviour of W,(@) to be non-singular along 
the axis of symmetry, and thus set A = B. Furthermore, we shall also 
exclude the region 6 = z in the present discussion. 

We shall next consider the solution /f,(@). Based on the previous 
discussion, we have shown that it is justifiable to write 


ie i 


> et ea —a)+4/3+f,0080. (4.10) 


This equation retains the same form as that studied by Yatzeev (5) and 
Squire (2) for the non-conducting fluid. 
In what follows, we shall discuss separately the two distinct cases 


(1) @€=6=é=0 


and (2) @=b = —2. 


The physical meaning of these solutions will be discussed in the following 
sections. 


Case (1) &@=b=é&=0 
In this case, the solution of f,(@) is simply 
2sin 6 
a A 4.11 
fo cos 6—k’ ges 

where k is a constant. The properties of this constant will be considered 
later. 

Now with the known zero-order solutions f,(@) and W,(@), we may 
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further find sie and - successively from (4.3) and (4.6): 


W(6) = ———, ig sin 6’ f (ih <a ree 


ay, — l)log(k—cos +7 —— = {(k—-eos @)log(k—cos @) +- 


ie eceide hati eet (4.12) 
where ¢, and ¢, are constants. 
If we require W,(@) = 0 (non-singular at 6 = 0) and (d?W/dé*),_, = 0 
(symmetry), we must have 
e+e = aa [(2k+-6)log(k—1)]+ low? for W,(0) = 


Sey am) _ 
€, = 2+— for (a. = 0, (4.13) 
whence é, = "e 2(k eet he N+ 
Hence, since W,(@) is known, the solution of f,(@) is obtainable from (4.6). 


This first-order differential equation is integrable and the solution may 
be written as 


f(0) = — eee, | Se wm, a0 
_ 2A®sin@ ff (k—cos6)* 
(cos 0—k)? sin 8 i 
+log{(k—cos 9)k—cos 9 1+-cos 6)! +C08 Oy 2/0e+1) 4 de, cos 6+-e,)] dé. 
(4.15) 


The integrals cannot be obtained in closed form and the final expression 
is lengthy: 


i sin? { — 
fi) = —cos 6)? | — oy 


+(1-+4-cos 6)log(1-+cos 6)+(k-+1)]+[—{eX(k—1)*+ 
+-fe,(k+1)(k—3)}log( 1—c08 #)—3(¢+¢)(1-+e*) ———— 


log 2—2. (4.14) 








[log(k—cos 0)*-1+- 


[(k—cos @)log(k—cos 0) + 


1 
1+ cos 6 
+-{he,(k-+1)(k—3)—fe,(k2+ 6k+ 1)}log(1 +08 4) |+ 
(k+-1)8(k+-3) log(k—cos @) 
4(k+-1) B,0080-+H(k-+1)5| 1+cos0 
1 log 1.008 9 cos 6 
meet 5 econ 
+-4(k-+ 1)[log(1+-cos 6)]*-+-3(k—1)* B, cos 0}, (4.16) 


a 


|+ 1(k-+-1)(k2—2k—7)B, cos 0+ 
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where 
8, cos @ = —log(k—cos @)log( 1 —cos 0) +-}{log(k—cos 6) |? +- 
(k—1) If ?--l J lf k-1 8 
+e Seta cc36) * \k—coso) 7" 


1+cos@  k—cos@ 
+k k+1 - 


1/k—cos6\2_ 1/k—cos 0\8 
+i) +4 Fy aria 


8, cos @ = log(k—cos 6)log 





2 
8, cos @ = log(1+-cos @)log en 


— }(1+-cos #)—}{4(1+-cos 4)}*?—}{4(1+-cos 6)}5+.... 


Therefore, we may, in principle, calculate as many terms as we desire. 
In a similar manner we obtain W,(@) and f,(@), etc. 
The numerical results of W,, fy, W,, and f, are tabulated on pp. 18-19. 


Case2. G=6 = —2 


In this case, (2.11) may be written in the form: 


dfy _ 


—— LB-+fecot a + Lore e— I 


sin?6 


Based on a similar consideration given by Squire (3), we obtain 





~ 1+cos6 cos $p log(cos 6+ 1)+-ysin 4p log(cos @+ 1) 


i sin 6 [2 ea ae ar an are 
0 eal ’ 


(4.18) 


where y is an arbitrary constant to be determined and p = 2é—1. This 
result agrees with Squire (3) though it is given in slightly different version. 
In the following discussions, we will also restrict our interest to the case 


fol47) = 9; (4.19) 


fo(0) = sin 6 (y—p)sin }p log(cos 6+-1) 
. (paoe8 | dos Tploc ose foie $p log(cos 6+-1) 


then 








sin 0 1+p? 


~ 14-008 6| p cot $p log(1+-cos #)—1 (19%) 


The constant p here depends upon the jet momentum (3). 
Using the known solution of f, and W,, we may determine the solutions 
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of W,(@) and ng in the forms 


; a ee 
mo= af ane f fees 





1 
cos 4p log(cos 6” + 1)—sin 4p log(cos 6” + 1) 


4 ° a 
ithe ocsran | aes = ‘Tens 8) | sin 6” | ee 


sin 6 





| aorae’, (4.21) 








sin 0” (1 —cos 0”) dé’”d6” ; 
* (1+c0s 0)" pcos 4p log(cos 6” + 1)—sin $p log(cos 6” +- mikes 
(4.22 
The boundary conditions on W, and f, may be satisfied by choosing proper 
integration constants. 

Although the closed form results of these integrals are difficult to obtain, 
we may proceed, in principle, to find the higher order solution as far as 
we like. For example, the second-order solution may be written in terms 
of the zero-order and first-order solutions in integral forms 

6 


e 
4 ] e , 7 , , 4 , Y Yr ” , 
W,(9) ain 6 | sin @ | [Wofi+Wfo—h Wo—fo W;)d0"d0’, 
fal@) = —sin Bel tead fs ~— [aw 24 We W—4f 3] 00’. 


5. Physical significance of the solutions 

At the beginning of this paper, we postulated the functional forms of 
the solutions H(r,@) and v(r, @) [expressions (2.5) and (2.6)] without giving 
any physical meaning for the choice. The analysis so far is thus mainly 
of mathematical interest. Now, the question naturally arises: What is 
the physical problem corresponding to the present solutions? To clarify 
this point, we examine more closely the previous results. 

In the first case evidently the zero-order solution gives Landau’s ‘un- 
bounded jet’ solution (1) with no magnetic interaction. The constant 
appearing in (4.11) may, therefore, be determined by the formula 


4 fil 


M = 16m tok + en) 8; (5.1) 


which was first discussed by Landau (1). - (5.1) M is the total impulse 
across a spherical control surface. This expression shows that for a weak 
jet (Mf — 0) the constant & will approach infinity and, on the other hand, 
k — 1 for a very strong jet (M — oo). The second case of zero-order solution 
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corresponds to Squire’s submerged jet diverging from a hole in a plane 
wall, 

Hence, an interpretation of the solutions which we discussed in the 
previous section is that they describe an unbounded jet of an electrically 
conducting fluid with an externally imposed magnetic field 

A’ (cos 0—1) A’ 4np\} 

aH, => —- ——-—— =_l A = ee A ’ 
RT ae ae Hro r ( be 
which may be considered as the field due to the distribution of ‘magnetic 
poles’ along the negative axis. In the previous section, we employed the 
solution scheme 
6 . . 

ne > W,(0) > f,(0) > Wy(8) > f_(8) > .... 
The physical meaning of such a procedure can be explained as follows. 
Conceptually, one may think of the fluid as primarily non-conducting. 
In other words, at the beginning we simply have Landau and Squire’s jet 
solution and we call that the ‘basic flow’ [ f,(@)|. Next, we postulate that 
an external field H, is imposed on the flow field. At this moment, since 
the fluid is non-conducting, no interaction between the flow field and 
magnetic field will occur. But, if the gas now becomes electrically-con- 


ducting, the ‘induced magnetic field’ [W,(@)| will appear because of the 
motion f,(9). However, this induced magnetic field in turn modifies the 
flow field to f,+af,(@). This new flow field will again distort the magnetic 
field so that W = W.4+aW,+o7W,. This mutual interaction will continue 
until an ‘equilibrium’ stage is reached. 


6. General remarks and conclusions 
a A’ A’ (cos@—1 
Since i, = — Hoy = rs ( — 9 LL 
the equation for the magnetic lines of force may be readily seen to be 
r(cos@—1) = constant. (6.1) 
The plot of these magnetic lines is given in Fig. 1. Now because H, is 
both irrotational and solenoidal, there exists a potential ¢ such that 
¢ = A’ logr(l+-cos@). 
The equipotential lines thus will be 
r(1+-cos 6) = constant. (6.2) 


For simplicity we shall postulate that the second- and higher-order 
solutions are negligible. In other words, the hydromagnetic interaction is 
significant only to first-order solutions of W and f. In case of an unbounded 
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jet, it is not difficult to show that the equations of the magnetic lines of 
force after interaction are 


r|(cos §—1) + 2a(k—1)log(cos 0—k) = [(cos @—k)log(cos @—k) — 


a 
k+1 


log(k—1)+ pare? — |} = constant. (6.3) 


k—1 
log ” o 


—(cos 6+ 1)log(1+-cos ay}+2a{ (1 an 


Joos 6— 


__ (k—1(k+3)4+2 
(k+1) 
The streamlines, on the other hand, can be found similarly by integrat- 


ing the following equation: 
rdé dr 


fotof; if Fy+aF,’ 
or simply rsin 6( f+ af,) = constant. (6.5) 
As f, and f, are known, the plot of such a family of streamlines is not 
difficult and is illustrated in Figs. 2, 3, and 4. Similar discussion may be 
extended to the case of wall jet. 

Another point of interest in the present problem is to investigate how 
the presence of a magnetic field shifts the isothermal lines of the non- 
conducting solution discussed by Squire (2). Then the temperature distri- 
bution may be found by solving the equation: 


eT eT e( 22 L Oh. oe - 
,OF Ol oe Nabeaatierecd wa dF 6 
or + 2 | ( o)+ wand salsin® =) iets 


where & is the thermal diffusivity and the dissipation terms are neglected. 
If we let 7 —T,, = [T(@)]/rt and use 


6 d 
ve = AM; ¢, = —*| foot 0+ $5, 








(6.4) 


we obtain 
hd ia 
ape Sites (6.7) 


where Pr = v/k = p/ck is the Prandtl number and c the specific heat. 
Since when » = 1, dI'/dy = 0, the solution I'(7) will be 


I(n) = R,e-Privia—nndn, (6.8) 


where R, is a constant determined by the strength of the heat source. 
If the fluid is non-conducting 
2sin 6 2(1— 7?) 
h= ak ANE Rane 4 
cos 6—k n—k 


~oosd—h’ * ” ond —E” 


+ JT. is the temperature at infinity. 
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In case when the gas becomes electrically conducting and hydro- 
magnetic interaction occurs, the temperature distribution will be 


T(r, 0)—T,,(r, 0) = B, e-Prfife+afi+..Wa—nildn 
r 


ie 1 

~ # (cos6—k)2Pr 
Thus, the new temperature distribution can be obtained by multiplying 
by a factor e*PrfSitafe+.0d0, 

Lastly, it is seen that the numerical computations indicate that the 
expansion procedure used for finding the previous solution gives fairly 
reasonable results. It is believed that the higher-order solutions are 
negligible in case « < 10-* for the present problem, although a rigorous 
proof of the convergence of such an expansion scheme, in general, appears 
impossible. 

In case (1), unbounded jet, 6 = 7 is excluded from the region of our 
interest. Therefore, the validity of the method of solution at this particular 
line 6 = 7m is immaterial. However, it is seen that the temperature distri- 
bution near the negative axis (9 = 7) is greatly affected by the presence of 
this singularity. From the numerical results, we can see that significant in- 
teraction appears in the region |@| > 42. This fact is reasonable since in this 
region the intersecting angles between the streamlines and magnetic lines 
become larger than that in the region || < 47 and also H, is stronger there. 

In case (2), wall jet, it is highly interesting to observe that maximum 
induced velocity occurs in the vicinity of 4°. This phenomenon is due to 
two facts: first, in the region of small @, the intersection angle between 
the #-component of velocity and magnetic field is large; secondly, it is 
seen from the calculations that f,(@) has a maximum somewhere between 
0° and 10°. 

All numerical results of f,(@) and f,(@) and W,(@) presented in Tables 1, 
2, 3, and 4 assume A = 1. 


e®Prf(fitafe+.d0 (6.10) 
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TABLE 1 
Case 1 (pipe jet), A = 1 





W, i 
Wo a pe rer [K = 3 € Ke=rijkK= (= 2)K=riiK= ror 








° ° ° ° ° ° ° 
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— 2°1446 o 8791 1°6475 7°3119 | 16°1753 | I°5571 | 7°2023 | 10°4750 
2°6558 06889 2°4067 | 10°1092 | 21°5906| 2°9571 | 11°7720 | 16°5800 
— 3°7321 0° 5088 3°8288 | 14°6558 | 30°3860 | 5°3467 | 21°1898 | 26°3105 
— 56713 0° 3006 6°4973 | 23°4971 | 47°5515 | 11°5343 | 42°7777 | 881510 



































TABLE 2 
Case 1 (pipe jet), A= 1 
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—0°6348 —0°9457 —0°9943 | 11-9148 | 46°7705 | 94°3567 

—0°6529 — 09495 —0°9947 13°0285 45°7114 85°5423 
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TABLE 3 
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TABLE 4 
Case 2 (wall jet), A = 
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THE FINAL STAGE OF DECAY OF A LOCALIZED 
DISTURBANCE IN A CONDUCTING FLUID 
IN A UNIFORM MAGNETIC FIELD 
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SUMMARY 


A general localized disturbance is given to an incompressible conducting fluid 
of infinite extent in a uniform magnetic field, and the final stage of decay is 
examined for the case in which a net linear momentum is imparted to the fluid. 
It is found that after a sufficiently long time the velocity field is that due to the 
superposition of two identical viscous vortex rings, of the type found by Phillips (1) 
for the corresponding problem in a non-conducting fluid, whose centres move in 
opposite directions along the uniform magnetic field with the Alfvén wave velocity. 
Some brief remarks are made about the decay of a turbulent wake behind a body 
moving through the fluid. 


1. Introduction 


THE motion generated when a general localized disturbance is applied to 


a (non-conducting) incompressible viscous fluid initially at rest has been 
considered by Phillips (1). It was supposed that the fluid was either of 
infinite extent or that the boundaries were so far away that they had 
negligible effect on the motion. It was further shown that, when the 
disturbance imparts a non-zero net linear momentum to the fluid, the 
motion in the final stages of the decay, when the non-linear terms in 
the equations of motion may be neglected, is a type of viscous vortex 
ring. In this the vortex lines are circles with centres on a line parallel 
to the net linear momentum vector and whose planes are perpendicular 
to this direction, the centre of the viscous vortex ring being at rest rela- 
tive to the fluid at infinity. 

Phillips applied his results to the final period of decay of the turbulent 
wake behind an axially symmetrical body moving through the fluid. In the 
present note we shall investigate the analogous situation for a conducting 
fluid in a uniform magnetic field. In the absence of such a field, there 
is no interaction in the final stage of the decay between the velocity and 
magnetic fields (supposing that a magnetic field was created initially along 
with the velocity field). The reason for this is that the coupling terms are 
quadratic, and the velocity field in the final stage is the same as if the 
fluid were non-conducting, although the motion in the initial stages of 

(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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decay will, of course, be different. If no magnetic field is created initially, 
it remains zero throughout the entire motion, that is the same as if the 
fluid were non-conducting. 

When a uniform magnetic field is imposed a different situation arises, 
since the velocity and magnetic fields are then coupled by the uniform 
field. It is now not obvious what the motion is in the final stage of decay. 
It will be found that there are then two viscous vortex rings, which are 
not at rest relative to the fluid at infinity. They move in opposite direc- 
tions along the uniform magnetic field with the Alfvén wave velocity 
corresponding to the strength of the uniform field. This result is made 
the basis for some remarks on the possible nature of the turbulent wake 
behind a body moving through the fluid. 


2. The equations of motion 


For a non-conducting fluid; it was pointed out by Phillips that the 
vorticity w(x,t) is exponentially small at infinity since a fluid element 
can only acquire vorticity by viscous diffusion, and that its Fourier 
transform 


Q(x, t) = (2n)-* [ w(x, the in.x dx (1) 


is therefore an analytic function of x and, in particular, can be expanded 
in a power series about x = 0. Let u(x,t) denote the velocity field and 
¢p(x,?) its Fourier transform; then Phillips showed that 


;(x, t) = (3.,—“23) 4+ 01%), (2) 


where 87%9M is the net linear momentum of the fluid and is constant. 
The right-hand side of (2) is not analytic at « = 0, although it is bounded. 
This is associated with the fact that in general the velocity field is O(r-*) 
at infinity (r = |x|) so that it is not absolutely integrable. 

In a conducting fluid, vorticity is also generated by the Lorentz force 
and it is now no longer true to say that the vorticity is exponentially small 
at infinity. Using e.m.u. the equations of motion in the magnetohydro- 
dynamic approximation are: 

div B = 0, curlH = 47j, curlE = -3; (3) 
j=o(E+uAB) B=uH; 
0. Vu = — yp + f AB; 
at p p 


divu = 0; 
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where B, H, E, and j are the magnetic induction, the magnetic intensity, 
the electric field, and the current; o, u, v, and p (all assumed constant) 
are the conductivity, permeability, kinematic viscosity, and density; and 
p is the pressure. From these equations H and w are given by 


= = eurl(u A H)+AV7H, (7) 
A 


where A = (47yc)~! and is called the magnetic diffusivity, and 


“Oo +-(u.V)eo = (oo. Vu-+rV%0-+_ curlj AB). (8) 
é p 


Equation (8) shows how the Lorentz force generates vorticity unless j \ B 
is the gradient of a scalar, which is not in general the case. 


3. The behaviour of the disturbance at infinity 

We assume that as r > 00 the disturbance vanishes so that the magnetic 
intensity becomes uniform and all other variables tend to zero. The 
irrotational part of the velocity field is derived from a scalar potential 
which satisfies Laplace’s equation, in virtue of (6); thus, in general, the 
velocity field will be O(r-*) at infinity. It follows from (7) that @H/ét is 
O(r~*), so that in general H—H, 


h (9) 
H, 


will be O(r-*). It now follows from (3) that j is O(r-5), and the contribution 
to éw/ét from the Lorentz force is O(r-*). Hence, in general, the vorticity is 
O(r-*) at infinity. (In the absence of a uniform magnetic field, an element 
of fluid only acquires a magnetic field by magnetic diffusion, so that the 
vorticity is exponentially small at infinity, as for the non-conducting case. 
It is essentially the coupling of the velocity field with the magnetic field 
at infinity which produces a vorticity field which decreases according to 
a power law.) 

The properties of the Fourier transforms of w, u, and h, near the origin 
of wave-number space, now follow from the behaviour of these variables 
at infinity. Since the first two integral moments of the vorticity field exist, 
the first two derivatives of its Fourier transform also exist, and hence 


€02,(0, t) 2Q,(0, t) 
j + rcs ey, ——— + 


CK; OK; OK}, 


0,(x, t) = 2,(0, )+K« o(x?), (10) 


where (see (1)) 
e02,(0, t) = LE sik M.. 


CK; 


Q,(0,t) = 0, 
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It then follows (again see (1)) that 


$,(x,t) = (5 —“32)m+ O(n), (11) 
so that the leading term in the expansion of the Fourier transform of the 
velocity field about the origin of wave-number space is, in fact, unaltered. 
(This is not the case for the higher order terms.) 


The Fourier transform of h(x, ?), (x,t) say, has the expansion 
T(x, t) = 1(0, t)-+0(1), (12) 
and the requirement that divh = 0 implies that «,T(0,t) = 0, ie. 


r(0,t) = 0. It follows that the Fourier transform of the current distribu- 
tion is o(x) near x = 0. 


4. The dynamics of the motion 


Choose axes such that H, is in the 1-direction, i.e. H, = (H,, 0,0), and 
denote the Alfvén wave velocity for this field strength by a, so that 


a® = pH?/4mp. 
It can be shown by straightforward but somewhat cumbersome tensor 
manipulation [for example, see Batchelor (2)| that the Fourier transforms 
of (5) and (7) can be written in the form: 


me +u«d, —a*ix, E == ine 3,5 8) Ay — a* B,.;), ( 13) 


a +AcP, -ik, b; = iK,(C;;- C5), (14) 


Cc 
where 


A ,;(x, t) = [ p;(x—x’, t)dh,(x’, t) dx’ | 


EB, (x,t) = geal (15) 


Cij(u, t) = ( $(x—m’, LT ;(x’, t) dx’ 


Note that A,,(x,t) is the Fourier transform of u,(x, t)u,(x,t), and similarly 
for the other expressions. The kinetic energy of the motion is 47°pA,,(0, t), 
and the magnetic energy associated with the disturbance is ?H? B,,(0, t). 

Since, in general, w,(x,t)u,(x,t) is O(r-*), it follows that A,; can be 
expanded as a power series in «x, as far as, and including, the terms which 
are O(«*). similarly B;, can be cxpinced as far as terms O(«‘), and C,,; as 
far as terms O(x*). Hence, substituting (11) and (12) into (13) and com- 
paring terms of order zero in «;, we find that @M,/ét = 0, i.e. M = constant. 
Thus the net linear momentum of the motion is a constant of the motion, 
as in the case of a non-conducting fluid. This result follows from the 
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conclusion that h is O(r-*) as r > 00, since the Maxwell stresses are O(r-*) 
and the net stress over the sphere at infinity vanishes. We can also deduce 
that the net angular momentum of the motion is constant, since the 
moment of the Maxwell and viscous stresses over the sphere at infinity 
also vanishes. [This result can also be deduced formally from (13).] 

It follows from equation (14) that the value of the magnetic spectrum 
near « = 0 depends on the history of the motion and is not a constant 
determined by the initial conditions. 


5. The final stage of decay 


The fact that the leading term in the expansion of the velocity spectrum 
function about x == 0 is a constant of the motion enables us to determine 
the asymptotic form of the velocity field as t-> oo. In the final stage of 
the decay of the motion, u and h will be small and their squares and 
products in the equations of motion may be neglected, so that the equa- 
tions can be linearized. The linearized forms of (13) and (14) are 

* +vx*4,—a*%ix, T, = 0, (16) 


i eT in $i = 0. (17) 


These equations have solutions proportional to e~-™, where 
(m—«*v)(m—x?A)+a%? = 0. 
The two roots of this quadratic are 
My = 4x2(A+v)+4b{xt(A—v)?—4a%P}}. (18) 
The roots are real and positive if jax,| < }«x?|A—v|; otherwise they are 
complex with positive real part equal to 4x?(A+-v), corresponding to damped 
periodic oscillations. 

Let ¢\°(x) and I'{°(%) denote the values of the velocity and magnetic 
spectra at some time ¢, within the final stage of decay, i.e. t, is a time for 
which the linearized equations are valid. The order of magnitude of ¢, 
depends upon the initial conditions and the history of the motion. The 
value of the velocity spectrum at subsequent times is 


(x, t) = e-tA+nt- °| [cosh 16-1) — "ih 40(t—t) 444) + 


+ — P\(*)sinh J01t—t)| , (19) 


where © = {«*(A—v)?9—4a%F}. (20) 
Now when (A-+-v)(t—t,) is large, the exponential term effectively damps 
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out the velocity spectrum except in the neighbourhood of « = 0, so that 
to a good approximation, and making use of (11) and (12),¢ we can write 
$/(2, t) = e-x*0+-mt-40 cosh 101-1) —0—) sinh 40(—t)(5—“3?) 4, 

(21) 
provided that the net linear momentum is not zero. The velocity spectrum 
given by (21) depends only on M and is independent of the history of the 
motion, so that the asymptotic form of the velocity spectrum, and hence 
of the velocity field, depends only on the initial value of the net linear 
momentum of the fluid. On the other hand it can be shown that the 
magnetic field in the final stage of decay depends on I'\°(x) and hence on 
the history of the motion, so that a ‘universal asymptotic form’ for the 
magnetic field does not appear to exist. 


6. The velocity field for the case A = v 


The expression (21) is still fairly complicated, and before examining the 
general case it is useful to consider the case A = v for which the analysis 
is greatly simplified. We then have © = 2iax,, and (21) becomes 


d; (x,t) = e-**vt- 6(8,—“) My, cos ax,(t—t,). 


The vorticity spectrum corresponding to (22) is 
Q(x, ) = te 454; Me" 'c08 ax, (t—t). (23) 
The solution of Phillips (1) for a non-conducting fluid corresponds to 
(22) and (23) with a = 0. 
The vorticity field, obtained by taking the Fourier transform of (23), is 
w(x,t) = —MAVO, (24) 


where 


(x,t) = |e te-#-0008 ax4(t—t9) dx 


7 \fexp| st lex] sag 


~ 2\v(t—ty) ~ 4y(t—to) | Letty) 


+exp|—| (25) 


ne nll 


4v(t—t,) |) 
When a = 0, this vorticity field was shown by Phillips to correspond 
with what he termed a ‘viscous vortex ring’. When a + 0, it is clear that 
the velocity field arises from the superposition of two identical viscous 
vortex rings. Each of these is the same as the viscous vortex ring which 

+ The second term in the curly bracket does not remain small compared with the first 


when a tends to zero, and is retained in order that the expression (21) be a uniformly valid 
approximation. 
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would appear in a non-conducting fluid when the net linear momentum 
of the disturbance is 4(87%p)M, except that now the viscous vortex rings 
are not at rest relative to the fluid at infinity but move with velocities +a, 
respectively, in the direction of the uniform magnetic field. The ‘radius’ 
of each ring increases like [v(t—t,)]* and the distance between them like 
2a(t—t,), which therefore becomes large compared with their radii. When 
this is so, each ring contributes independently to the kinetic energy, which 
then decays like [v(t—t,)|-!. It is to be noted that these moving viscous 
vortex rings do not carry the fluid along with them but are essentially a 
wave-like motion in which the vortex lines move through the fluid with 
the Alfvén wave velocity. 

One point worth making is that the present result is unexpected, since 
it appears that the sole effect of the uniform magnetic field is to split into 
two the viscous vortex ring which would appear in a non-conducting fluid. 
Intuitively, one would expect the magnetic field to damp out the velocity 
components across the field whilst the velocity component along the field 
was relatively unaffected, but this is not the case. 


7. The case A ~ v 
In the general case, the asymptotic vorticity field is again given by 
(24), but now 


° 2 
@(x,t) = | e-t*A+t-t0)+ {cosh 10(t—t,)—G sinh ,0(¢ -ty)} dx, (26) 
where © is given by (20). An evaluation of this integral in closed form 
does not appear possible, but it can be evaluated approximately by the 
method of steepest descents for large values of t—t). 

The main contribution to (26) comes from those values of x for which 


* 


in.x—}(t—t){(A+v)x2?LO}, 


regarded as a function of x, is stationary. Omitting the analysis, which 
is straightforward but cumbersome, it can then be shown that if 


a*(t—ty) ~ 
“A—y 


me (27) 


then the method of steepest descents gives a value for ® identical with 
(25) except that 4(A-+-v) replaces v. That is, when (27) is satisfied, the final 
stage of the motion consists of two viscous vortex rings, each of which is 
associated with a momentum 4(87%)M, moving with velocities +a in the 
direction of the uniform magnetic field and decaying as if they were in 
a non-conducting fluid of kinematic viscosity }(A+v). This result appears 





FINAL STAGE OF DECAY OF A LOCALIZED DISTURBANCE 27 
to be true whatever the value of A/y, although this ratio will affect the 
value of f,. 

2(¢__ 
If, on the other hand, : (t—to) <1 (28) 


jA——P | 


for some time within the final stage of decay, then the velocity field is 
the same as if the fluid were non-conducting. That is, a single viscous 
vortex ring occurs which is at rest relative to the fluid at infinity. How- 
ever, after a sufficiently large time (27) must always be satisfied, so that 
the single ring must eventually split in two, although by the time this 
happens the kinetic energy of the motion may be extremely small. 


8. Remarks on a turbulent wake 

We conclude with some speculations on the nature of the turbulent wake 
behind a body moving sufficiently rapidly through a conducting fluid in 
a uniform magnetic field. For a non-conducting fluid, it is found in practice 
that the wake is turbulent whenever the Reynolds number is sufficiently 
large. The turbulent wake may be thought of as arising from a distur- 
bance to the fluid caused by the body. Sufficiently far behind the body, 
where the Reynolds number of the turbulent wake is small, the motion 
can be regarded as due to a superposition of viscous vortex rings. Now 
for a conducting fluid in the presence of a uniform magnetic field, it is 
not absolutely certain that a turbulent wake exists since the magnetic 
field may prevent the transition to turbulence. However, this is most 
unlikely, especially if the Reynolds number is sufficiently large and the 
magnetic field is not too intense. We may then expect a turbulent wake 
to oceur which will decay as the distance away from the body increases. 
The present analysis then shows that in the final stages of decay the wake 
will split into two. Relative to axes fixed in the body, at a sufficiently 
long distance from the body, the wake will consist of two parts each of 
which is similar in structure to the wake when the fluid is non-conducting. 
It is inclined to the velocity of the fluid at infinity, the inclinations 
depending upon the orientation of the magnetic field and the ratio of the 
Alfvén speed to the speed of the fluid at infinity. If this ratio is greater 
than unity, and a turbulent wake still occurs, then one branch of the wake 
may extend upstream of the body. 

Although the present analysis only shows that the wake in the final 
stage of decay splits in two, it is probable that it splits in two before the 
final stage is attained. The reason for believing this is that turbulence 
is perhaps best recognized by the existence of a randomly fluctuating 
vorticity field. Now one can regard an Alfvén wave in an incompressible 
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fluid as a vorticity wave, so it does not seem altogether unreasonable to 
suppose that some of the vorticity associated with the turbulence is 
propagated along the lines of force with the Alfvén wave velocity. How- 
ever, it does not seem possible at present to say more than this, and it is 
hoped that future work may elucidate this problem. 
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SUMMARY 


In their linearized theory of the axisymmetric supersonic free jet Taunt and Ward 
(1) encountered a certain function. However, they did not tabulate it or consider 
its detailed properties. 

In this note the properties of a generalized form of their function are considered 
and a short table of the original function is given for the range of argument 0 to 8, 
in intervals of 0-1, to five decimal places. 


1. Introduction 


TAUNT AND Warp (1) treated the linearized theory of the supersonic 
axisymmetric free jet by an integral equation method and found the 
boundary shape of the jet to depend on a certain function, denoted by 
F(z) by them, whose Laplace transform was J,(p)/pJ,(p), where the symbol 
I, denotes the modified Bessel function of the first kind of order n. Ward 
(2) (in a more accessible source) has given the theory in an entirely opera- 
tional form. 

In both the above references only a graph of the function was given; 
a complete investigation of the properties of this function or a table of 
its values was, however, not given. 

The present author has had occasion to compute a short table of F(z) 
and to investigate the properties of a generalization of this function (3). 
These may be of some general interest. 


2. The functions F(z) 


A generalization of the function of Taunt and Ward (1) is considered 
heret and is denoted by the symbol F,(z). The functions F,(z) are defined 
by the operational relationships 


_ LP) 
* pl, (p)’ 


F, (2) = (1 ) 
+ This work forms part of a thesis for which the degree of PhD. of the University of 
London was awarded. 
{ This function can be shown to arise in the analysis for a jet with non-uniform boundary 
pressure. 


(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961) 
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using the symbolism of Van der Pol and Bremmer (4). The function F, 
corresponds to the function F mentioned above and will be denoted hence- 
forth by F. 


3. The properties of the functions F, (z) 
By the inversion theorem of Laplace transforms, we have from (1), 


, l LAG 
L 
where # is an infinite path to the right of all the singularities of the 
integrand. 

It follows at once as a standard result for one-sided transforms that 
F(z) = 0 for z< 0. For z > 0 the contour may be completed by an 
infinitely large circular are to the left of # which is approached by a 
sequence of finite arcs not passing through a singularity of the integrand. 
As |I'.(p)/p*I,(p)| is of order |p|-? for large p, the integral along this arc 
vanishes in the limit and so F,,(z) is given by the sum of all the residues of 
the integrand. 

If jm is the mth positive zero of J,(z), the ordinary Bessel function of 
order n, the integrand has simple poles at p = +4%),,,,- Also j,,9 is zero 
except for n = 0, but there is, nevertheless, still a simple pole at the 
origin when n = 0, produced by the factor p?. 

It is then found that for z > 0 


F,(z) = nm +4nz*—2 > COB Jn? 
2(n+1) ~* = Fam 


For large values of m the following holds [see Watson (5)|: 


. l 
jn ~ Ya2n+4m—1)+0(= (4) 


Therefore the series on the right-hand side of (3) is uniformly and absolutely 
convergent, since the moduli of the individual terms are ultimately less 


than or equal to S11 = 
and fm TEE oe 
at 7? 6 
m=1 


mn? 
Also, by inversion of the asymptotic series for the integrand, the 
following expansion for z > 0 in powers of z is obtained: 
F,(z) = z—}2?+ O(2). (5) 
For n = 0 this has been extended to give 
F(z) = z—0-252? — 0-0208333z3 — 0-0052083z4 — 0-001627 625 — 0-00056422* — 


—0-000207 927 — 0-000079825 — 0-00003162°—0-0000121z!°—.... (6) 
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The series (5) and (6) certainly converge for |z| < 2; z = 2 being the 
position of the first of the singularities of F’,(z) as will be shown. As 
F,(0) = 0 from (5), (3) then yields 


ont) See, 
~ 4(n+1)’ 


= (7) 

Jam 
an interesting result which, besides being useful in computing the F,, 
also indicates that (from (3)) 


Se 
dnz? < F(z) < gnzt+ —- (8) 


i" 
Thus, F,(z) is continuous everywhere and is bounded for finite z. There- 
fore, provided that the resulting series is uniformly convergent, (3) may 
be differentiated to give, for z > 0, 
F,(z) = nz+2 > SID J nm = (9) 
m=1l1 Jn m 
F’,(z) has singularities and discontinuities which are now investigated 
by a method used by Ward (6) in a similar context. This also shows the 
regions in which the series (9) is valid. Using (4) 
SID J nm pa sin(m— }+ }n)xz 10 i ; 
pea mr m* 
Now the sum of terms O(1/m?) is uniformly and absolutely convergent as 
before and since for convergence only the ultimate behaviour of the series 
for large m is of consequence it follows that 





(10) 





Re (SIN J»m% sin(m—}+ }n)r2) 
22. oo "Uee e (11) 


mr j 


nym 


is a uniformly and absolutely convergent series. Thus any singularities or 
discontinuities of F),(z) are contained within the function 


“ S sin(m—}+4n)xz 


mr 





m=1 


2{.. > cos maz —* sin mz) 
= — )— Z ——_—— + —h)haz — . (13 
aaa $)4r 12d — cos| (n — }) 47 1D ae (12) 
The two series on the right of (12) are readily summed? and it is then found 
that F(z) has the same singularities and discontinuities as the expression 
. 2 i 
os sin{ (n— })42z}log[ 2(1—cos mz) ]+ — coal (n -- i)drftan(; ons} ; 


(13) 





x 
‘ : C08 mnz ++ sin maz 
+ Since —log(l1—e'*?) = ° 
m 


m=1 
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This is subject to one qualification. The expression (13) is algebraically 

valid for z < 0 and is symmetrical about z = 0, whilst F’,(z) is only given 

by (9) for z > 0 (apart from the points of singularity and discontinuity), 

being zero for z < 0. Thus only half of the discontinuity of +-2 at z = 0 

which is possessed by (13) is realized in F’,(z) itself—a conclusion entirely 

in agreement with the derivative obtained at z = 0 from the series (5). 
Thus it is found that F',(z) has singularities of the form 


2 (—)™+*log z—48s— 2} (14) 
7 


at z=48+2 (s = 0,1,2....), 


and discontinuities of the form 
(—).2 


at zen de (e = 1,32,3.,...), 


with, in addition, a discontinuity of +1 at the origin. 


4. The computation of F(z) 
Setting n = 0 in (3) yields 


« 


F(z) = 3-2 > Sein’ (16) 
m=~1 Jom 

This equation, together with the power series (6), forms the basis of the 

computation. The asymptotic formula (4) for j,,,, shows that 


COS Jom? Cos(m—f)nz 


17 
ris > (17) 


and so, generally, the series does not converge rapidly enough for direct 
calculation (to five decimal place accuracy) to be profitable. However, 
for certain specific values of z, such as z = 2, the cosine term leads to 
more rapid convergence and (16) may be used directly. 

For some other values of z the series may be transformed to a more 
rapidly convergent form. The most important instances are when z has 
odd integral values or when it is a multiple of 4, }, or any power of }. 
It is then found that the signs of successive terms rapidly become those 
of the ultimate form (17). If now, without changing the order of the terms, 
they can be taken either singly or in small groups (which are separately 
summed) to form an alternating series of slowly decreasing numbers, then 
the well-known Euler transformation can be used to obtain a more rapidly 
convergent series. For the smaller sub-intervals obtained in this manner, 
however, even this method becomes too lengthy because of the number 
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of terms of the basic series needed to form each of the terms of the 
alternating series. 

This method was used to compute F(z) at intervals of } up to z = 2 and 
the series (6) (which agrees perfectly with such values in the region z < 1-1, 
where the number of terms shown in (6) is adequate for calculation to five 
decimal places) was used to obtain values at intervals of 0-1 up to z = 1-1. 
It may be noted that the useful range of the series (6) can be somewhat 
extended by forming from it the power series for 

F(z)—(2/m)(z—2) log (1 — 32). 
The values of F(z) at points symmetrically placed about z = 2 with respect 
to the foregoing were then found, using the expression 


cos Zo, = ©CO8 Jom? 





F(2+2)+F(2—z) = 1—4 S 


m=1 Jom 


(18) 


which follows from (16) and converges rapidly enough for practical pur- 
poses. 

The framework so obtained was subtabulated to an interval of 0-1 in 
the region 1-1 << z < 2-9 by use of Langrangian interpolation with the 
interpolating function 

F(2+2)— F(2—z)—(4/m)zlogz (19) 
and with (18) also. 

The region 4 < z < 8 was readily covered using the formula 


cos COS 4) om m ©O8 Jo. m2 


F(4+2)+ F(4—z) = 1—4 "i 


Jam 


(1+-cos 4}6,m)©O8 Jom 2 


2[1— F(z)|—4 (20) 


m=1 Jom 
which also follows from (16) and is rapidly convergent in the last form. 
Expressions similar to (18) and (20), which converge rapidly, may be 
easily found to extend the existing values beyond z = 8 if required. 


5. The table of F(z) 


The table given below can be interpolated directly except from about 
z = 1-7 to 2:3 and from z = 5-7 to 6-3. In these regions the interpolating 
functions 


F(z)—= (2—2)log|1 —42| and F(z)+ = (2—6)log|1—4e| 


can be used, except very close to the singular points when (19) together 
with (18) must be used near z = 2; and similar expressions must be used 
near z = 6. It should also be noted that as F’(z) is discontinuous at z = 4, 
separate difference tables should be formed above and below this point. 

D 
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F(z) F(z) F(z) F(z) 


0"00000 , 029668 ° 087468 . 076979 
009748 . 022273 . 076641 . 0°83065 
018982 4 018399 , 066804 p 085412 
027689 , 016766 057942 . 085441 
035851 “s 016808 , 050025 6 0°83798 
0°43451 . 18083 : 043038 , o81001 
0° 50466 ; 020484 , 0° 36963 077204 
0°56874 ; 0°23793 ‘ 031795 ré 0°7 2666 
0°62646 . 027831 "27524 0°67596 
0°67748 3° 0° 32490 . 024158 ° 0°62130 
O7 2143 





075781 3 0° 37696 5° 021709 
078620 3 0°43375 . 020184 
0°80589 3° 049466 5: 019608 
081547 3 0°55909 ; 020081 
0°81320 3 0°62654 5 0°21730 
079855 0769648 . 0°24554 
076784 3 076850 P 028850 
71619 38 084219 5° 0°35033 
063212 og1719 . O°44131 
046108 “ 099315 rl 0°61277 


0° 56362 
0750385 
0°44265 
o° 38081 
0°31890 
025756 
0°19728 
013854 
0°08170 
0°02714 
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Fie. 1. The function F(z). 
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SUMMARY 
A solution is obtained for the steady flow of a visco-elastic liquid past a sphere 
at small Reynolds number. The drag is found to be less than that given by Stokes’s 
formula for a Newtonian liquid of the same viscosity. 


1. Introduction 

In 1958 Oldroyd (1) considered non-Newtonian liquids whose rheologicai 

equations of state relating the stress tensor S;, and the rate of strain tensort 
Ey, = HOpg t+ Uy) 

take the form Sy = Px—P oy, (1) 


P**+-),(d/dt)P*+-y, B® Pi+-v, E*P,g* 
= 2n | E*+-A,(d/dt)E“+-v, EX E,g*}|. (2) 
Here U; denotes the velocity vector, g,;, the metric tensor, and P an iso- 
tropic pressure; 7) is a constant having the dimensions of viscosity and 
Ax, Ass Mo» Vp» Mg aT constants having the dimension of time. The derivative 


denoted by d/dt is the ‘convected derivative’ (Oldroyd (2)). For any 
tensor B,,, we have 


s Bik — M Be. UB. Wi, Boks We, Bim Bi Bmk_ Bk Bim. 
where ia = Ws — Uy). 

Oldroyd found that in theory liquids of the class defined by (1) and (2) 
exhibit non-Newtonian flow properties which have been observed in visco- 
elastic liquids. Such liquids have a variable apparent viscosity in simple 
shearing, decreasing with increasing rate of strain from a value 7, at low 
rates of strain to a value 7, (< 7m) at high rates. They exhibit the 
Weissenberg climbing effect and have a distribution of normal stresses 

+ Covariant suffixes are written below, contravariant above, and the usual convention 


of summation over values 1, 2, 3 applies to repeated small Latin suffixes. A suffix ¢ following 
a comma indicates covariant differentiation with respect to the space variable 2’. 


[Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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associated with an extra tension along the streamlines in many types of 
simple flow, having anisotropic state of stress normal to the streamlines. 
Also for flow at small shear stresses, such liquids are characterized by 
three constants, a coefficient of viscosity no, a relaxation time A,, and a 
retardation time A,. 

If the above properties are to be obtained for all rates of shear, Oldroyd 
found that the following restrictions on the six constants in equation (2) 


are necessary : e >'s, >be, (3) 
where 01 = Ay pot (Ar— Fo) 

and Og = Ag uot (Ay—FHo)r2 

and (Ay —Bi49)(Ay ¥—Aav) > 0. 


In this paper a solution of equations (1), (2) and the familiar equations 
of motion and continuity for steady, incompressible flow of a liquid of 


density . 
ensity p Si, = pU*Ut,, (4) 
EFi=0 (5) 
is obtained for uniform flow past a sphere, the solution being a perturba- 
tion of that due to Stokes. A power series expansion in terms of the 
dimensionless group UA,/a is found as far as the second power in that 


parameter, where a is the radius of the sphere and U is the velocity of the 
uniform stream. The solution is valid if 


[FS 


<i, 

No @ 
ie. if a® <A, n/p and U < a/A,. Thus there is an upper limit to the radius 
of the sphere and the stream velocity. When , = 10 P, A, = 10~' sec, 
and p = lg/e.c., a? <1 em? and U < 10a cm/sec. For larger n. the 
restrictions are less severe. 

The pattern of streamlines obtained differs very little from that found 
by Stokes for a Newtonian liquid. Upstream the streamlines are slightly 
displaced towards the sphere and downstream they are displaced away 
from the sphere. 

The drag on the sphere has no linear term in UA,/a but contains a 
quadratic term in this factor. Furthermore, it can be shown with the 
aid of (3) that the coefficient of this term is negative. Hence for a visco- 
elastic liquid in which the relationship between stress and rate of strain 
is given by equations (1) and (2) the drag on the sphere is less than it 
would be for a Newtonian liquid of viscosity np. 

Values for the relaxation time A, and the retardation time A, for a 
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particular liquid can be obtained by experiments with a cylindrical visco- 
meter (Oldroyd (3)). However, the author knows of no estimates for the 
other three constants jy», v,, and v,. If sufficiently accurate experiments 
could be performed with ball-bearings in a visco-elastic liquid the results 
obtained would provide a measure of ¢,—o,. This certainly would not be 
sufficient to determine each of jo, v,, and v, but would be a way of esti- 
mating their magnitude. 

Recently Tanner has carried out some experiments with ball-bearings 
in a visco-elastic liquid and gives a report on his findings in an appendix 
to this paper. There is good agreement between theory and experiment. 


2. The equations of motion 

Spherical polar coordinates (R, 0, 4) are chosen with origin at the centre 
of the sphere and @ = 0 in the upstream direction. Since there is axial 
symmetry the solution is independent of the coordinate ¢. All tensor 
quantities are expressed in terms of their physical components referred 
to spherical polar coordinates. 

The equations are reduced to non-dimensional form by the substitutions 


R = ar, U, = Ua,, Up = Uugy, 
= 1o(U/a)pP, Poe = nolUla)pep, Pog = no(U/2)P gg, 
Pe = 1(U/a)pe, P= n(U/a)p, 
E,, = (U/a)e,,; Egg = (U/a)egg, Eug = (Ulaegs, 
Eg = (Ula)e, Wo = (U/a)w,g. 
Also it is convenient to write 
As=A,, pe =A, y, = £A,, ve = mA,, 


where ¢, {, €, 7 are dimensionless constants. Equations (2), (4), (5) become 


Prt] ra. FPn + “Hl; Pn 2a) + Plt —¢9)— 2 P| + 


+ Le re( Pre t+ Poe + Pop) + £7 (Cor Pre +€00 Poot € 44 Pog + 2€r9 Pra) 


= 2 eg bert tnt = —2eq) + 2¢,9(w,9—e,9)-- 2e2, _ 


+ nr(e3,+-¢55+-e34+ 2%), 
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Pee + “| mez = Peet: “(5 7p Poet 2p) — 2p,9(W9+€r9)— 256 Pa + 


+ Cre g6( P+ Peet+ P44) +ET (Cre Pret €00 Poo + ph Pop + 2€r0 Pro) 
= 2lean berm ui, <tgp us “a(S eo0+ 2) ~ Zeta +e) —2ey| + 
+ ane tely techy + 20)], 
Pest t| ex = Pat 8 Pes 244 Pas] + 
+ bre 44( Pre t+ Poot Po) +E7 (Cre Pre +00 Poe +€4¢ Pog + 2€r0 Pro) 
- 2legg-ter| u Seas + awe] + 
+a toy ty t 2c, (8) 
pat*| us 5p Prot = M(t Pr Pos) + Bol —€19) — Pr Wp + lg) — 


— Palen +¢e)| 4 Cres(Pre+ Peo +P4e) 


= 2eeter| w= —€rg+ “5 COTE —¢e) + altro —€0)— 


—tnltoa te) —etlen te) |} (9) 


1 @é 1 op 
naw fh oleies ate —(2p,.—pyg— cot #) —— 
Pret = appt = ( Prr—P00—P 46+ Pro ) Or 


or 


— Pata < Poot ~{(Poo—Pyg oot 0+3p,4}—— 2 


= Rel ue te? — utp 

eu, a+ ; Fy hn 

1 @ 1 

—~ — (r*u,) +-—_—. 

73 or! 4 rsin 6 60 
UDA, 


where T= —, 
a 


e (up sin 0) = 0, 


Equation (12) shows that a stream function % can be introduced such 


that 
fee 1 
“ ™ — tend 00’ “o™ 5nd or (13) 


The boundary conditions are that the velocity of the liquid is zero 
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at the surface of the sphere and U is in the direction of 6 = 7 at large 
distances from the sphere. Hence 
u-=U=0 atr=1 

and u, = —cos8, Up =sin@ atr = o. (14) 

In order to obtain a solution for the above equations it is assumed that 
the Reynolds number Re is small and therefore the right-hand sides of 
equations (10), (11) are neglected. This implies that Re <r, i.e. 

a® < Ai No 
ee 

and hence the size of sphere is limited by this relationship. 


3. Solution of equations 

A solution is obtained in the form of a power series in 7 for quantities 
P, Up, Ug, Derr POO Phd» Pro> Ps Ever 00s Ebd> Sng» Wyg; 1.€. 

A = A®M+47AD+472A@+ ...,, (15) 

where A stands for y, u,, Wg, Ppp, P99 Pd PrO> Ps Crrr €00> Cbg» rg» Wyg in turn. 

If the expressions (15) are substituted in equations (6)—-(13) and the 
boundary conditions (14), and terms independent of + are equated, the 
following equations are obtained: 


ha = 2Qe(9, - = 2h, PY = 2B. pip = 29, (16) 


0) 
(0) (0) (0) __ (0) __ (0) (0)) _ op 
4! ran +- “(2p Poe Pog t cot Opy) oe ty ae a7) 
7 
= 2 + e 4s {9 — pif} )eot 0+-3p9) — 1 op 
apt 90 r 0 


¢ (sin Ou) = == 0, (18) 


2 1 
p2y() 4. _ 
a or )+ rsin 6 20 


1 oy (0) l 
—-- ,F7 Uu = —_— -—» 19 
resin? 00 © 0 rsin@ or (19) 
uo =u =0 atr=1, 
ui) — —cos8, up = sin@ at r= oo. (20) 
Using (19), (16), and (17), @ép/ér and (1/r)(ép/20) can be expressed in 
terms of Y and its derivatives. ga p™ an equation for Y is 
obtained, viz. ae sin@ @ af * is ” 
dr?" 7? A0\sin @ 5 Gs 7 
The solution of equation (21) for boundary conditions (20) is 


YO = 4(r?— r+ 1/2r)sin®d (22) 





pear kemfee tne <Not aati fap Ae eas See anlar 
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as found by Stokes (4). Iv is now possible to obtain all the terms inde- 


pendent of 7 in expressions (15). 


If expressions (15) are substituted in equations (6)-(13) and boundary 


conditions (14) and coefficients of r equated then 


Ata a 


ans 8.8. Ay. 
7 (\-5-a+3-za)}+ 


4-2 (2% (15) + 


\ 
ies 2. 5 “) 


" 
Qe) — x Grote. ote th clara amb 
™ ‘)) r eA Ss 


oe 


2 1\2 
(a 


2 
(1-2-3 
Ane. 


1 1 
atest aal|+ 


2 


4% —#)/5%8 “(1—3)'+ 


269+ ae 
i rr pS 


r 


(6 9 16 18 +3) 


e 


. ani 2 ny! 
or 


or 


1 @ 1 
= aii — oe eee 


crm nya t Proc 


r sin 3 iil? as 


oy uD — iieae 1 


~ -sind 66” 9 rsin 


“=u =0 atr=landr =o. 


ayy 


rsin@ @r’ 


sin’ 
r8 


mn, 
re | 


sin*6) 
— | 





: ag Pre +. t (2p? —piy —p\) + pPeot 6), 


—PY)+3p(9}, 











42 F. M. LESLIE 
From equations (23), (24), and (25) 


op é 2 2 2 
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sin @ ay) erase. r® 














A. 
r=| r=2 6=0 
Fic. 1. The streamlines y) = constant ; i.e. the secondary flow due to non- 
Newtonian terms to the first order. 


The solution of equation (29) for the boundary conditions (28) is 
3 
yd) — H(1—e)sin*# cos 0(1—*) ‘ (30) 
. 


From (30) it is possible to obtain all the coefficients of + in expressions 
(15). 
The streamlines ¥ = constant are shown in Fig. 1 for the flow in the 
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first quadrant only since there is axial symmetry and the flow in the 
second quadrant is the mirror image of that in the first quadrant. The 
direction of flow is towards the line @ = 47 in both first and second 
quadrants. 

By exactly the same procedure as for ¥ and ¥ the following equation 
is found for ¥, 


@ |, sing sin@ @ e\|* yo 
or? * 7? BO Te 


_ 2) 2 {sin 27 160 236 75 148 
™ oe (22-5 r + r* r5 
5541, 46, 50 56)) 
era ie a ae | 
4 (1 —¢) (in? 27 440 424 75 = 138) | 
Fe(1 os wel. eee or ae ear 
4sin*6/), 5 aa 86 50 54 Pe 
re a ee 
sin 40 243 240 120 75 42 60 
(oo-28 rm ra i 


+f(l—e oe 


9 
yp_202, 78, 42 100_ 12 58)) | 
r r2 r? jf 


HEA 8 att I—atate 


where 8 = (n—£)(1—8£)—{(1—e). The solution of (31) for the boundary 
conditions 
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A PCS r ro 
From the inequality (3), 
ATS +-(1—BL)E} > ATSe+(1—FL)y}, 
i.e. (1—$L)(n—€)—C(1—e) = B < 0. 


Thus £ = (o,—o,)/A2 and is always negative. The streamlines 
J = constant are shown in Fig. 2 for « = } and 8 = 0 and in Fig. 3 for 
« = 4, 8 = —4 for flow in the first quadrant only, there being symmetry 
about @ = $7. Here, however, the flow is always in the downstream 
direction. The streamlines for « = 4, 8 = —1 were also obtained but 
show little change from those for « = }, 8B = —}. It is found that at 
r = 3, 0 = }r, with « = }, 


, 03429 logr 0-1609 00-1905 0-0154 | 


ZO: J: J® = 1 : 0-029 : 0-014—0-4068. 


Thus the convergence of this series for the stream function would appear 
to be satisfactory provided 8 < 1 and + < 1. Since the ratio of ¥ to 
YW is large the streamlines % = constant are little different from those 
for ¥ = constant if r < 1 (ie. UA,/a < 1). As a result in Fig. 4 the 
streamline % = 0-373 is sketched alongside that for Y — 0-373 but the 
difference between the two is exaggerated for clarity. . 
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rel r=2 6-0 


Fic. 2. The streamlines *) = constant fore = 4,8 = 0; i.e. the secondary 
flow associated with the square of the parameter UA,/a. 
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ri r=2 6=0 
Fic. 3. The streamlines i) = constant fore = }, 8 = —}; i.e. the secondary 
flow associated with the square of the parameter UA,/a. 
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Fig. 4. A sketch of the streamlines = 0-373 (with « = }, 8 = —}) and 

) = 0-373 comparing the solution for a visco-elastic liquid with that for 

a Newtonian liquid. represents (i.e. Newtonian). 

represents ~% (i.e. non-Newtonian). The difference between the two is 
exaggerated. 


4. The drag on the sphere 
The drag on the sphere is given by 


D = —2na* | (S,,),-18in 0 cos @ d0-+2na* { (S,9)--1 8in6 d8, 
0 0 


nT 


Eee 0 = | (P,@)p-1 8in*6 dé +- [ (@—Pn)--r8in 8.0089 dé. 
2any Ua 
i) 0 


If J, J, / are known, it is possible to calculate the quantities required 
to determine the drag. It is found that 
D 
—__. — 3—[0-016(1—«)(3—«) —0-6188}r?. 
ag 3—[ 6(1—e)(3—«)—0-6188]r 
This is equivalent to 


v3 


D = 6g Va ——"2°— [0-016(A, —A,)(3A, Aq) —0-618(,—o,)}. (32) 


Hence to first order in + the drag is unchanged from its value for a purely 
viscous liquid but to the second order the drag is decreased, since 8 < 0 
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and « < 1 (ie. A, > A,). Unfortunately no values are available for 8 so 
that the magnitude of the decrease cannot be estimated. 
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APPENDIX 
By R. I. TANNER 


The general conclusions derived in the above paper have been tested by dropping 
steel balls of various sizes (# in. to 74 in. diameter) through a solution of polyiso- 
butylene in carbon tetrachloride. The value of yn) was approximately 90 P, and the 
diameter of the fall tube containing the liquid was 40 mm. The ball was timed over 
a distance of 150 mm, split into two sections. In many cases observations over 
both halves of the fall were taken to provide a check on the homogeneity of the 
liquid. For each size of ball the mean of a number of fall times was taken, care 
being taken to drop the various sizes in a sequence which minimized the effect of 
temperature variations. Mean fall times varied from 70-2 to 21-5 sec. 








7 ‘ 

(U/eP (sec*) 
Fic. 5. Preliminary experimental results for steel balls in polyisobutylene 
solution. 


It was calculated that the maximum Reynolds number Re was about 0-006. 
This would provide about a 0-1 per cent. deviation from Stokes’s law in the New- 
tonian case, and was therefore neglected, since it was completely masked by the 
wall effect. According to B.S8.8. (5) the increasec drag (D) due to the wall proximity 
is given by the following formula (for a Newtonian liquid) 

D = Dj (1—2-104a/h+ O(a*/h*)}], (Al) 
where D, is the Stokes drag, a the ball radius, and h the fall tube radius. The 
cubic term was negligible. 

In the present case the method adopted was to calculate the apparent viscosity 
(m) for each ball size using the measured fall time and expression (A 1) and extra- 


polate to a = 0 to find my. (D,—D)/D, is then equal to (y)>—7)/n» and has been 
plotted in Fig. 5 against (U/a)?*. 
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According to the theory presented above (equation (32)) 
D,—D 
where K = 0-005(A, —A,)(3A,—A,)+0°206(0,—o,), provided UA,/a < 1. Considera- 
tion of Fig. 5 shows this to be the case for (U/a)* < 2-5. The expected value of A, 
for this liquid lies between 0-3 and 1 sec so that the agreement is very satisfactory. 
It is possible that some of the deviation for large values of (U/a)* is due to wall 
effect which would alter the ‘elastic’ part of the solution, or to breakdown of the 
polymer under increasing rates of shear. These preliminary results are being ex- 
tended in the hope that they may provide a useful method of characterizing visco- 
elastic fluids. 
The author wishes to thank Dr. G. Allen of the Chemistry Department, Man- 
chester University for providing the materials for the polymer solution. 


= K(U/a)*, 
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POLYNOMIAL REPRESENTATION OF SOME SHIP 
SECTION AREA CURVES AND THE CALCULATION 
OF THE ASSOCIATED WAVE RESISTANCE 
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SUMMARY 


In this paper methods are discussed for fitting polynomials to empirical data for 
application to the theory of wave resistance of ships. The study arose from the 
discrepancies observed between resistances computed by use of various Lagrange 
interpolation polynomials. The limitations of this approach are explained here, and 
an alternative method of curve fitting by least squares using Chebyshev poly- 
nomials is described which provides coherent estimates of the resistances. A com- 
parison is also made with a method based on the Fourier analysis of the data. 


1. Introduction 


Ir is customary for some purposes in naval architecture to specify a ship 
in terms of ‘section area curves’. These = ___— = 


curves represent the area y of vertical sec- 
tions of a ship’s hull perpendicular to a line y} 
joining the bow and stern, as a function d 


of the distance x of the section from an 
origin on the bow-stern axis (see Fig. 1). In 
the present context we are concerned mainly 
with the rate of change of y with respect 
to x, and we therefore restrict ourselves to 
curves from which the central portions of 
constant area have been removed. The ee 

scale of coordinates is chosen so that the sisataba 

bow and stern are the points (+ 1,0), (—1,0) respectively. We refer to 
the range 0 < x < 1 as the forebody, and the range —1 < z < 0 as the 
afterbody. 


Hogben (1) has shown that, according to Michell’s theory, the wave 
resistance is given by 

















Ship's Hull 





R= [ (M?+ N*)cosh*u exp|_% coshu| du, (1) 


i 


0 
(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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_ {xg L cosh u 
2v? 


(2’) 


Here V is the speed, Z the length, and Z the mean draught of the ship, and 
g is the acceleration due to gravity. We see that the resistance correspond- 
ing to any prescribed section area curve y = y(x) may be regarded as a 
function of two parameters V/vL, L/Z. 
Equivalently, we have 
1 


+ 
M = a L oth | y coe *) dx, (3) 
~j4 


+1 
s L xg L cosh , 
N= i 40osh w | ysin( | dx, (3’) 


= 
since y= Oata = +1. 

In general the section area curve is specified by numerical values 
accurate to three or four decimals at the discrete points x = —1 (0-1) 1. 
The concern of this paper is to determine, according to mathematical 
criteria, a reliable and convenient method of calculating the value of R 
to be associated with such imprecise data. The problem is important since 
it is not possible to measure wave resistance directly, but only the com- 
bined effects of wave resistance and friction. 

One approach to the evaluation of R would be to use the numerical 
values of y directly in some appropriate quadrature formula for (3) and 
‘’). However, the inaccuracies which arise in methods of this kind are 
not readily assessed, and the possibilities of the approach have not been 
explored. 

The procedure investigated here is to approximate the data by an 
analytical expression in the form of a finite sum 


y(z) = 3 a, F(a) (4) 


of convenient functions F., and to substitute this in the expression for 
R. For example, if F(x) = 2’, then from (2) and (2’) 


V2 


r=1 


+1 
n e n 
M => ra, | e-tsin( OO) ae = 34,75, say, (5) 
» pis 
1 


arena te Memes ete, ata tl cenbirit ie brane! ean 
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= xg L cosh u , 
and = ae pms “¥eo9( HE pa ) dx = 3 a,rC,_;, say. (5’) 


Hence on squaring these results and substituting in equation (1), we find 


that 
= 3 3 4,0,hy (6) 


r=18= 
where 


is) 


L,=re | (S,_,S,-1+0C,_, C,_,)eosh*u exp|—fzcoshtu| du. (7) 
/ 


The computation of R may thus be separated into the evaluation of 
the J,,, which depend only on the parameters V/vL, L/Z, and the deter- 
mination of the coefficients a,, which depend on the given section area 
curves. The latter stage is the principal subject of this paper. 

Three cases are considered : (i) F(x) = 2", and y(x) determined by colloca- 
tion with the data, (ii) F(x) = cos(rcos~'x), and y(x) determined by least- 
squares approximation, (iii) F(z) = sinjrmz, and y(x) determined by 
collocation. 


0 


2. Choice of polynomial representation 

A simple method of obtaining a representation of the curve y = y(z) 
is to caleulate the polynomial of least degree which has collocation with 
y(x) at a set of discrete points: this is the Lagrange interpolation poly- 
nomial. Table 1 gives a selection of values of the resistance computed by 


TABLE 1 


Values of R corresponding to Lagrange polynomials 


Model r Model 2 Model 3 


30 7o 30 ° 30 7o 





0030 = 0060 oo016 }§ 0034 o088 860206 
0026 ©8605 0024 0053 o0g8 §=—: 0" 220 


059 «O22 0064 oOl12 0535 0992 
0054 OvITO 0069 «= O11 33 o-601 1104 


1'2 Ll, 0250 060208 0324 0479 0467 0-980 
L, 0244 0°383 0324 0496 0522 = r'089 

+ Here and in Table 4, the factor 1-689 arises from the measurement of V in knots and 
of L and Z in feet. 





this method for three section area curves illustrated in Fig. 2. L, is the 
value of RF calculated from the Lagrange polynomials which fit the fore- 
and after-bodies separately at the points |z| = 0, 0-4, 0-8, 1-0, and L, is 
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the value of R corresponding to collocation with the curves at |x| = 0, 
0-2, 0-4, 0-6, 0-8, 0-9, 1-0. There is a discrepancy between L, and L, which 
varies from less than 1 per cent to 30 per cent and is generally greater than 
10 per cent. 





as On 


“+ \ Model 1 
-N 


Modei 2°. \ 
awa, 





i 





Fie. 2. Section area curves examined. 


These large disagreements suggest that the available data have not been 
used properly in approximating the section area curves. However, it is 
undesirable merely to increase the number of ordinates which are fitted 
since this would entail an increase in the degree of the Lagrange poly- 
nomial, and it is well known (2) that the exact fitting of many data values 
introduces another type of inaccuracy—namely a violent oscillation about 
the curve which represents the true function. 

A superior procedure is the approximate fitting of all the data by a 
polynomial of low degree. For example, the method of Least Squares 
yields more satisfactory results and in addition provides criteria on the 
adequacy of the fit. The theory of least-squares polynomial approximation 
is well established (see, for example, (2)), and its practical application using 
an automatic computer has been considered by Forsythe (3). The follow- 
ing modification of the method is due to Clenshaw (4). 

Let y(2,), y(%e),--., y(%,) be the data values. Then the polynomial 
approximation y,(x) of degree n which minimizes the residual sum of 
squares, 


Fh = {Yn()—Y(@)}P +{Yn(%2)—y(%2) Ph +--+ £{Yn(2m)—Yl(Zm)}?, (8) 

is evaluated in the form 
y, (x) = dy+d, T,(x)+-d, T,(x)+...+d,, T,,(x). (9) 
T(x) denotes cos(rcos-'x), the Chebyshev polynomial in x of degree r, 
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and the coefficients d, are obtained without solving sets of linear algebraic 
equations. In general, the set 7\(x) (r = 0, 1,...,m) is not orthogonal with 
respect to the points 2,, 2»,..., 2,,, 80 that the effect of passing from y,(z) 
to y,,.,(@) is to change all the existing coefficients dy, d,, d,,..., d,. The 
magnitude of these changes and of the new coefficient d,,, is a useful 
guide in assessing the ‘goodness of fit’. This is an alternative criterion 
to the variance-ratio test (described, for example, by Quenouille (5)) which 
tests the hypothesis that d, ,, is zero by examining the magnitude of the 
variance ratio F,,,, where 
,  &_,—8 
F, = 2-1, 


on 


(10) 


and o? is the residual mean square 52/(m—n—1). 
We now describe the application of this method to the present problem. 


3. Results of least-squares fitting with Chebyshev polynomials 
The section area curves were represented by values of y at the 21 points 
x 1(0-1)1 for models 1 and 2, and at 18 points for model 3 (for which 
readings were not available at x = —0-3, +0-5). Least-squares poly- 
nomials of all degrees up to the seventeenth were calculated in 10 minutes 
for each model using a Deuce computer, a machine having a multiplication 


time of 2 msec. In addition to indicating the ‘correct’ fit, the results 
identify data values with outstandingly large residuals over successive 
polynomials. These values impede the smooth fitting of the aggregate of 
data and constitute ‘unfairness’; they are therefore checked for errors. 

The Chebyshev coefficients d, of successive fittings for model 2 are given 
in Table 2 together with the largest positive residual «,, the largest 
negative residual e_, and the corresponding values of x. It will be noticed 
that the even coefficients are unchanged in passing from an even to an 
odd fit, and vice versa. This is a consequence of the orthogonality of the 
auxiliary polynomials used in constructing the fit when the data points 
are symmetrically distributed. 

On examining the trend of the coefficients in the Chebyshev series and 
the magnitudes of the largest residuals, we might infer that the fifth degree 
polynomial would be an adequate representation. However, the statistical 
evidence presented in Table 3 suggests otherwise. The residual mean 
square o,, diminishes perceptibly with increasing n until n = 10 and then 
remains approximately constant, and the variance ratios F, to Fig are 
significant at the 5 per cent level (6), while F,, and onwards are not. This 
indicates that the tenth degree polynomial is required. 

If the residuals for the two polynomials are compared, half the fifth 
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TABLE 2 


Chebyshev coefficients of least-squares polynomials 


Degree of 
polynomial, n 





Chebyshev coeffi- 
cients 


+0°4010 +0°4012 +0O°4012 + 0°4003 +o'4011 
—0°0269 — 00269 -O°0251 — 00263 —0°0263 
—0°4940 — 04938 —0°4938 —0°4954 —0°4940 
+0°0420 +-0°0420 +-0°0432 -+-0°0422 +0°0422 
+0°0935 +0°0931 +-0°0931 +0°0923 +0°0932 
O'O103 —O°0103 OO1l2 oro1r16 —oOo116 
-0°0026 — 00026 —0°00290 —O°O0I7 


VEER BES Ee 


oR the 


— 00080 -o" —0°0074 —0°0074 
+ 0°0054 +0°0046 
+ 0°0031 +0°0031 
dig — 00032 


Largest residuals 
and data points 
€, 0°0077 


+o9 


0°0096 
+07 





TABLE 3 


Statistical evidence relating to least-squares polynomials 





Degree of Residual mean Variance ratio, 
polynomial, n square, O, F 





1°6 
ogo 
oo! 
0°42 
o19 
o12 
O-044 
0046 
0040 
0046 

















degree residuals are about ten times greater than the known uncertainties 
in the data, and groups of consecutive residuals have the same sign. The 
size of the tenth degree residuals is close to that of the data errors, and 
like signs occur consecutively at only four pairs of points. The choice of 
polynomial depends, therefore, on whether smoothness or closeness of fit 
is considered the more important. 

This uncertainty is also present in calculations for the other models. 
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Accordingly the wave resistance was calculated for the polynomials sug- 
gested by the two criteria and also for those of intermediate degree. 


4. Calculation of the wave resistance 
On replacing y in equations (3), (3’) by the right-hand side of equation 
(9) and then squaring and substituting in equation (1), we obtain with 
the least-squares polynomial of degree n 
Rw > > 4d, Ky (11) 
r= P = 
where, in place of the J, used in equation (6), 
K (S, S,+C,C, . (_9% . 
n=\(sp a) fe i jeosh wexp| — 73 cosh u) du, 


0 


_— - | ‘epinf Dosh dx, (13) 


ae f T,{x)00s| eee (13’) 


9V2 | 


We may note in passing that the substitution of equation (9) in equations 
(2) and (2’) would lead to a different expression for R, since y,(x) does 
not vanish at z = +1. This alternative expression for R is likely to be 
less accurate because essentially it involves the numerical differentiation 
of an approximate formula. 

For computational purposes we write 


e~"(1+-2/b)[C, C, +S, 8] dt, (14) 


+1 
S, = | T(x)sin{xa /(1+-t?/b)}¥da, (15) 


=-j 


C, = [ T,(x)cos{xa,/(1-+-t2/b)} dex. (15’) 

1 
Infinite integrals of this form are evaluated conveniently by special use 
of the trapezoidal rule (7), C. and S, being calculated for each value of ¢ 
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by the following recurrence formulae : 


: ee 4n., 2 sindA 


potana <Being an > We ak a: ’ 
4 cosa 
2n—3 A 
2sinA 48; 


eas ieee Ge 


Con -2 + 


a ae 
= = (n>2), (16) 


Ing seat 
Cy 2cosA 
ictal Sali Weg 
where A =: a,/(1+-¢?/b). 

Evaluation of these integrals constitutes the bulk of the work of 
calculating the resistances. However, the K,, may be tabulated against 
V/VL, L/Z once for all, and the resistance corresponding to any section 
area curve is then obtained as the sum of weighted products of poly- 
nomial coefficients. For a polynomial of degree eleven, for example, this 
takes about two minutes using the Deuce. 

The resistances R,, calculated by this method are listed in Table 4. 


TABLE 4 
Values of R corresponding to least-squares polynomials 
y : Model r Model 2 Model 3 


ccummemeuhinn 30 7° 3° 7° 30 7° 
1-68ovL 





og 0°039 

0°035 0100 

oo41 0-100 

or044 0°097 

0044 0097 

0043 0°096 
0-096 


0°063 

0062 o504 I'Ios 

0063 0594 I°I04 
8 0-064 0594 I°00Q7 
9 0-064 o594 I°097 
10 0-064 °o°594 1og! 
11 o595 x1°og! 





There is general agreement among the answers given by the different 
degrees, much more so than between the earlier estimates L, and L,. 
The discrepancy between each R,, for the extreme values of n nowhere 
exceeds 10 per cent, and for the larger value of the parameter V/vZ is 
less than 5 per cent. The values of R, decrease with increasing n for 
models 1 and 3, but the reverse is true for model 2. 

We may also note the similarity between the R,, corresponding to each 
even » and the succeeding odd n; this reflects the predominantly even 
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form of the section area curves. An exception occurs for n = 6 and 7 for 
model 2, R, departing considerably from the general trend. This illus- 
trates the danger of relying on a single fit of predetermined degree, but 
otherwise appears merely to be a consequence of the particular resistance 
terms corresponding to d, being relatively large, and not offset by higher- 
order contributions as in the subsequent polynomials. 

In brief we may say that although it is desirable to examine the values 
of R,, corresponding to a range of n, fairly accurate values may be 
obtained for the parameters considered by using the smallest n satisfying 
one of the criteria given above. 


5. Application of Fourier analysis 

Another method of curve fitting with evenly spaced data is described 
by Lanczos (8). In this approach we calculate the coefficients B, of the 
finite Fourier sine series 


g(x) = B, sin + B, sin 72+ Bysin > + ..., (17) 


which fits the data exactly at corresponding locations in the transformed 
range (0,2). The coefficients corresponding to an analytic function should 
decrease asymptotically as n~* in these circumstances, whereas the coeffi- 
cients corresponding to superimposed noise will not converge. Accordingly 
we expect the calculated coefficients to diminish rather steeply to relatively 
smal! values and thereafter to oscillate. For practical purposes, therefore, 
(17) is truncated where this oscillation commences. 

This method was applied to models 1 and 2, but although the coefficients 
decreased in magnitude they did not ultimately oscillate in sign. For 
model 1, B, to B,, were all positive, while B,, and B,, were negative. 
For model 2, B, to B,, were negative, and B,, to By, were positive. In 
view of this unsatisfactory result the method was not pursued. 

The weakness of this approach is that the data can be made periodic in 
the original variable z only at the price of a discontinuity of the second 
derivative at the ends of the range. Expansion in Chebyshev polynomials 
is also a Fourier expansion, but the independent variable is now cos~! x. 
The original function is thus transformed into a periodic function without 
discontinuity, and the convergence of the expansion is thereby greatly 
improved. 


6. Conclusions 
The method of least-squares fitting using Chebyshev polynomials has 


definite advantages over both Lagrangian interpolation polynomials and 
Fourier analysis. Firstly, we obtain a polynomial of low degree which 
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approximates closely all the data. Secondly, the convergence of the 
Chebyshev series is superior to that of the Fourier series. 

We note also that the rate of convergence of the Chebyshev series 
indicates the ‘goodness of fit’ of the approximation independently of the 
size of the residual mean square and the variance ratio, the usual statistical 
criteria. In the examples studied an examination of the Chebyshev coeffti- 
cients indicated a satisfactory fit of lower degree than did the statistical 
evidence. The wave resistance corresponding to both sets of polynomials 
were in adequate agreement, however. 
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DISPLACEMENT EQUATION OF ELASTIC 
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SUMMARY 
This paper is concerned with infinitesimal transverse displacements of a homo- 
geneous isotropic elastic material occupying the region of space —h <z<h, The 
essential feature of the analysis is the replacement of functions f(z, y,z) by their 


h 
moments (1/2h) j filx, y, z)z" dz (n = 0, 1,...). Simplifying assumptions are avoided 
-h 


until a stage in the analysis at which the errors introduced by such approximations 
can be estimated. The important question of the order of accuracy attainable with a 
fourth-order differential equation is considered. 


1. Introduction 


Tue difficulty of dealing with the exact equations of infinitesimal elasticity 
theory in the case of plate-like material has led to the treatment of plate 
problems by approximate methods. Mathematically the theory falls into 
two distinct sections: generalized plane stress, dealing with displacements 
which are essentially parallel to the plane surfaces, and transverse de- 
flexion theory, dealing with displacements which are essentially in the 
perpendicular direction. 

The theory of generalized plane stress is approximate to the extent that 
it does not serve to determine the actual stresses, etc., but only their mean 
values through the plate. In the case of transverse displacements, Green 
(1) has given a general solution for plates with stress-free surfaces, but 
application of the method to particular problems requires a considerable 
amount of labour. Also, it is not obvious what order of approximation is 
involved if the Fourier series or power series in this solution are truncated 
after a finite number of terms. Green and Zerna (2) assume that the 
stresses f*!, ¢? are quadratic functions of the distance from the mid-plane 
and derive an equation for ‘weighted’ displacements. The stresses across 
rectangular elements through the plate are then replaced by their statically 
equivalent force-resultants and couples at the mid-plane. Since this 
notion of equivalence is complete only for the application of dynamical 
principles, an approximation is here involved which leaves the order of 
accuracy of the final equations in doubt. 


(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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Other accounts of plate theory, making assumptions concerning stresses 
and/or displacements are given by Love (3), Coker and Filon (4), Musk- 
helishvili (5), Stevenson (6), (7), inter alia. Whatever the initial assump- 
tions, these theories are approximations in the sense that only mean values 
and first moment means (see, e.g., Stevenson (7)] are determined. This 
criticism does not, of course, apply to Green’s general solution. The most 
frequently made initial assumptions are those attributed to Navier or 
Kirchhoff [see e.g. Stevenson (7), Novozhilov (8),Green and Zerna (2), (9)]. 
In nearly all cases of this type of stress it is concluded that for normal 
loading with pressure p the transverse deflexion w, of the mid-plane 
satisfies the equation V‘w, = kp, where k is a constant. Love (3) realized 
that this equation is not accurate and obtained a correcting term. Green 
and Zerna (2) have obtained a similar type of equation for the ‘weighted’ 
displacements, containing a correcting term of like form. The significance 
of these terms is, however, not clear. 

An additional assumption that is sometimes made is that the displace- 


x 
ment vector u can be expanded in the form } u,(z, y)z", where z is the 
n=0 


distance from the mid-plane containing axes Oxy. Only the first few 
terms of this series are assumed to be important. There is some doubt as 
to the validity of this expansion and of its convergence. Further, as was 
emphasized to the author by Professor Green, the coefficients in the series 
will, in general, depend on the thickness 2h. In fact, the evidence in 


Green's paper (1) is thatu = Y¥ u,(z,y)(z/h)", where the u,,’s are inde- 
} k ~~. n . 


pendent of h. Thus, truncating this series does not readily enable the 
accuracy of the solution to be estimated in terms of powers of h. 

The present treatment is based on the fact that a function f(z) of a real 
variable z, satisfying sufficient conditions regarding continuity, is com- 


h 
pletely defined by its moments (1/2h) { f(z)z" dz (n = 0,1,2,...). Earlier 
and 


investigations have taken into account only the values n = 0, 1, a limita- 
tion which necessitates further assumptions in order to obtain a soluble 
set of equations. Here exact equations are retained until a late stage in 
the analysis. For a required order of approximation in terms of h, the 

0 
form of the equation for the mean value of the transverse displacement w 

0 
is seen to depend on the order of magnitude of w relative to h, information 
which can be obtained from the crudest theory. In general the order of 
0 

the differential equation for w depends on the accuracy required and, from 
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a practical point of view the most important question considered here is 
the accuracy attainable with fourth-order differential equations. No re- 
striction is placed on the conditions of stress over the plane surfaces and 
general body-force terms are retained throughout. 


2. Fundamental equations 
The fundamental equations of infinitesimal homogeneous isotropic 
elasticity theory may be written in the form 
t).+F° == 0, Fr = p( Fr—f"), (2.1) 
g* = 2ye* + Adg*", 2e* = g”"u®,+g*u" », (2.2) 
6= ge = gu, = w",,. (2.3) 
Thus te” = pig u®,+g*u")}+Ag* ul. (2.4) 
In the above equations ¢” and e* are the symmetric stress and strain 
tensors respectively, F’ is the body force per unit mass, f* is the accelera- 
tion vector, p the density, u’ the displacement, A and y are the Lamé 
elastic constants and , indicates covariant differentiation using the 
Christoffel bracket connexion formed from the metric tensor g,,. Elastic 
material is assumed to occupy the region —h <z <A referred to a 
rectangular cartesian coordinate frame Oxyz. For generality, curvilinear 
coordinates z', 2* are chosen in the mid-plane of the material. Thus 
xt = g#\(x,y) and 2* = x*%(z,y). The cartesian z-coordinate is retained 
throughout and is denoted by 2°. The line element ds is given by the 
2 ase ds? = g,,dx'da* = g,gdx* dx? +-(dx3)?. 2.5) 
Throughout this paper Latin indices take the values 1, 2, 3 and Greek 
indices the values 1, 2. Thus g,, has the following components: 
Jng(", x*) + @. Jas = Ian = 0. Iss = l. 
Also the determinant 
Ire! = \Gagl = 9 # 0. 


The contravariant metric tensor g* is given by 


g"° = 92/9. 9* = —guig, 9 = —Gulg, 9* = Gulg, 
whilst g” has the additional terms 
g@ = g™ = 0, g®@ = }. 


For coordinate transformations from z', z*, 23 = z to #, #, P =z 


ox 
—— #9 


ae f , 3 ar 
exr8 ; |) , 


és 
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and there are similar equations for the inverse transformation. Thus the 
transformation 
(2.11) 

for a contravariant vector gives 

ox* 

a* = —, uv, “2 = x3, 2.12 

x8 ne 

and there are similar equations for a covariant vector. In fact 
Up = Gry? OF ty =Gagw, ty =u. 
Similarly, for a double contravariant tensor field t* the relation 
Ox? Cat 

may be written in the form 


er exh orF* 
fap — té = 3 
dat ox0 os oat e 


2.15) 
iS — 233, we — 2 yoy 
ox 
Thus for the transformations considered, t** is a two-dimensional tensor, 
3 and #8 are vectors and f* is an invariant. The index 3 indicates 
invariance. 


3. The Christoffel symbols and the covariant derivative 
From the definition 


ial = $9°"(Ginjk + Gamj—Gikm)> (3.1) 


in which the comma denotes partial differentiation with respect to space 
coordinates, and equations (2.6), (2.9), it is easily seen that 
f3\_fa\_fo) 9 
Ne ath Fe an 0 (3.2) 
GR 3) (8 
i.e. all Christoffel symbols containing a 3 are zero. Returning to the trans- 
formation (2.10), the relation 


(3.3) 


between the components of the covariant derivative of a contravariant 
vector is equivalent to the equations 


— OX™ Ox" 


Uo = uf, ‘ 
B= Gat a 


“?,. — 2 
Ug = U'}s, 
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Similar equations hold for the covariant derivative of a covariant vector. 
Thus the index 3 again indicates invariance. For the covariant derivative 
of a double contravariant tensor field the equation 
My eat oe (3.5) 
Ox? Oxt OX 
yields the following equations: 
—p _ 0% axF dat ez aa? 
et tT a a 
(8, az8 ean Ox ag 
~ @xk ogy 
pas~ _ OF es a 6 = al 75 _ 433 
ak” [3° a Fe, /3 is J 
(3.6) 
As before, the index 3 indicates invariance. This is further emphasized 
by the form taken by the covariant derivatives. For example, the equation 





iia (3.7) 


gives the equations 
ux u* 41% uy; u%_ = a*,: ue = u® 
B po ly | ’ /3 ,3? ~ ix 


Similar equations hold for covariant vectors. Also, from the equation 


as ui, = u®,. (3.8) 


ow = Pot 


a (3.9) 


q p) q p) 


it follows that 
— 708 fo \ 988. ( B \ pad 
7 by 5 y| 
(3B, = BQ; tet, ae get 
“ jm ” , (3.10) 


= /33,.- 


3? #8, wi 





= 5 =, 
It is important to note that, since the Riemann—Christoffel tensor of the 
space considered is zero, the order in which repeated covariant differentia- 
tion is performed is immaterial. Also, since all Christoffel symbols con- 
taining a 3 are zero, parallel transfer of tensors along any line 


x! = const., xz? = const. (3.11) 


necessitates no change in the tensor components. Thus integration of 
tensors along these lines is permissible and further, since the components 
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of the metric tensor are independent of z°, all tensors, etc., derived from 
the metric tensor alone behave as constants in such integration. 


4. The moment equations 


n 
The nth moment F(z", 2?) of a function f(z',z*,x*) is defined by the 
equation 


e A 
Jict, 2%) = = [ fe, 2%, 29)22)" des in = 0,1,...). (4.1) 
<. 


The factor (2h)-' is introduced for dimensional convenience. From (4.1) it 
follows that 


n 


af » = of £-5 
- = —-f, —_ = hr-) — 9 4.2 
Cx iat x3 Rh af 9 


where 2f = f(z! 2?, +h) +f(z!, 2*, —h), (4.3) 
(+) 


the plus or minus sign being chosen according as n is odd or even. By 
definition the nth moment of the covariant derivative of a contravariant 
vector is given by h 


[ ul ,(x*)" da’, (4.4) 


From (4.4), (3.8), 
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Thus for the elements with Greek indices the operations of taking the nth 
moment and of covariant differentiation are commutative. Returning to 
the transformation considered in (2.10) and the vector transformation 


law of (2.11), the operation of taking the nth moment yields the following 
a 


h 
= l ox = age 
ae a’ (z3)" — « tg n i” = — yx. ‘ 
= ii’ (Z*)" d# Th f aa” (x3)"da*, or @ hs (4.7) 
—h —h 
Thus the operation of taking the nth moment has not altered the vector 
sith Similariy, from (2.14), 


Ox" ox 
r8( n —_ tP9(23)" dz, 
af? (#3)"dz#? = ae (x3) or 


so ae the tensor property of ¢”* is preserved. The result is quite general 
and holds when covariant differentiation is involved. For example, 


n 
= aFat— ; OR" 02 da —— 
Vy = —— 


—— uP ‘ . 
dx? axe t= top bat og 


(4.9) 


5. The nth moment of the fundamental equations 
It will be convenient to distinguish those terms which contain the index 
3 from those which do not. Equations (2.1), (2.3), and (2.4) will therefore 
be written in the form: 
Atty t+FP = 0, M% ty +F? = 0 

mF = pigruh,,+gPrur,}+dg%, 3 = uw", +u%s }. (5.1) 

$9 == Qpurg+Ad, 9 = HF = plu) +977U%,,} 
Taking the nth moments of the above equations and using the results of 
the previous sections the following equations may be obtained: 

n— n—1 n 

(., +h*- 10 —n 4+F8 = 0, 
=. aot 
(P = wl grat, + 9Pat, dl aaa u¥i, 


n—1 n 


) 
La (A+2u)[h- a? c @ | +Au",, 
+) 





* n- n—1L 
5) + hn ys3_n +7 = = 0, 
(B,) n n <*) n 5.3) 
pa = pa = ulin tus—n a +g, ° 
{+) 


Equations (5.2), (5.3), with n = 0, 1,..., replace (5.1) for it is known that, 
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subject to certain conditions regarding continuity, functions are com- 
pletely defined by their moments. The boundary conditions over the 
plane faces have been absorbed into the equations. The sets of equations 
(A,,), (B,) fall into two independent groups, viz. (A,,), (B.,.,) (r = 0, 1,...,) 
and (A,,,,), (B;,) (r = 0,1,...). The former correspond to stress systems 
of the type known as generalized plane stress and the latter to transverse 
deflexion of the plate. Here only the bending theory is considered but 
a similar investigation may be carried out for the generalized plane stress 
equations. 


6. The transverse defiexion solution 


Consider the two groups of equations (A,,,,), (B,,), viz. 

js ® etl 

PF + h*t8—(2r+1)88+ FB = 0, (6.1) 

as = ast tr 

8 = plo uP ), +9 we | +999 wr), thu? —(2r-+1)u'], 
(6.2) 

wit ew 

a — (A+ 2) (hu —(2r-+1)u8| +A ues (6.3) 

(+) 





\ 


2r 
(33, 4 por—1ys_ LF = (6.4) 


(-) 


(By) )e x ee 
| == f* = » giit<pteal u* +g*"u5,,}. (6.5) 


2r+1 
Eliminating f? from (6.1), (6.2) yields the following equation: 
ar+1 ar +1 ar 
(A+)g®* a ra ge” uP teat Ag’®| h?ru®—(2r+1)ui} +2198 
(+) (+) 
S SS 
—(2r+1)#?+ 48 =0. (6.6) 
— 
The value r = 0 presents no difficulty in (6.4), (6.5) as the moment ?¢* , 
2r 
though undefined in this case, is multiplied by r. Eliminating t** between 
the two last equations and also between (6.5), (6.6) gives the following 
two equations: 
ta = = ~ 
2r fF = pj h® tu — 2r us gut.) +h®-U34F°, (6.7) 
(-) (=) 
ev et a 
(A+ )g?* u* wat OO uf feat AgePh* us), +h2e8 + FB ae 
(+) (+) 
2r-1 ar 


—(2r+ ph 1uP + 2r(2r+ Le we —(2r+1)(A+p)g*us,, = 0. (6.8) 
(-) 
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Taking the covariant derivative and writing 
7? jap == A, 
(6.8) becomes 
2r+1 2r 


Al (A+2u) wB pt AiPru®—(2r+ 1+ my} + 


2r+1 2r—1 
+h¥ Pig + FB pp— (20+ 1)phY rahi g + 2r(2r + 1) w? ig = 0. (6.10) 


For r = 0 (6.7) becomes 
0 0 
p| Wet hud] ++AF* = 0. (6.11) 
{—) (—) 


Henceforth the following abbreviations will be used: 
a=A+2p, b=3A+4u, g=Atp, g=2A+3u. (6.12) 
The invariant uF ig will now be eliminated from (6.7), (6.10) by means of 
(-) 


(6.11) to give the equations 
2r r—1 ar 2r—1 1 ar 
ar = = pA| us — _prri\— our a ub gt F*— n#s, (6.13) 
ae A J wt 
Ala uw ptAhtud—(2r-+ 1)gud-+(2r-+ 1h?) + 2r(2r-+ 1) uP ig+ 
gia U 
+ FB pth O84 (2r+1)h*-1| 34+hF3\ = 0. (6.14) 
(+) Diy | 2r+1 
If r is replaced by r+ 1 in (6.13) the following equation is obtained for BS: 
ar+l 2r+2 ar+1 ar+2 9 
(r+1) = pal ut us _ pens —: Qu(r-+1) wg +F® —h*2F, (6.15) 
2r+1 
where r takes the values 0, 1,.... The moment f** may now be eliminated 
between (6.3), (6.15). Thus 
2r+1 ar+2 2 0 


2(r+1)g we s= = pA| we —h+ *u3) + 
ar ar+2 


+2(r+1)a{(2r+1 jus— hw) +F* _prrsnZs, (6.16) 
2r+1 
The moment wu can be eliminated between (6.14), (6.16). Thus 
2r+1 2r—1 
2(r-+ Val eetat FF wt Br 2r+ 1) uP jg) +2(r-+1)(2r-+1)gh®-1 x 
0 or+% 2 ar+2 2 
x [unser +094 nF) +aA|pAw +F% — ne ane + #4 
2r 


+2(r+1)(2r+ l)gAu—2(r+ 1)ubh* Aus = 0. (6.17) 
(+) 
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When r = 0, (6.17) becomes 


& 1 2 
2g{ +h] F + Hg + FP] + 2hbp{ w— us} ie 


2 2 0 0 
+ahA|pAu®+F%—h[pAud+F]| =0. (6.18) 
Replacing r by r+-1 in (6.17) the following equation is obtained: 


2r+3 2r+1 


2(r-+- 2) {O30 y+ FE ig + (r+ 1)(2r-+3)e wW ig) —2(r+ 2)ubh®r*Au + 


2r+4 ar+4 ar+2 


+aA(uA w + F% — ions Fe} +2 2r-+3)ugA uw + 


$2(7-+2)(2r-+-3)gh*r+ wh Aut +08 A = 0 (6.19) 
(-) 
at 
where r = 0, 1,.... Eliminating w* jg from (6.19) by means of (6.16) gives 
the following equation: 
2r+3 


2(r-+2)g besa Big + FP ig + (2r-+3)he 1) —2(r-+2)ubh*r Au? + 


2r+4 ar+4 
+aA{yA w + F > — — Weel ad + #')) ig 
ne 2 


+2(r-+2)(2r+3)y{ F +(2r+2)af (2+ 1) — hu?) + 
ar+2 
+4(r+2)(2r+3)yad uw UW +2(r+-2)(2r-+-3)Ahe+| Au + FA) =0. (6.20) 


After operating with A on (6.20) the invariant u® may be eliminated by 
(+) 


means of (6.18) to give the following equation involving only moments of 
the vector u*, which is henceforth denoted by w: 


abu 2 “ 
a a asl rs Iho —(r + 2) w }+ 


4 ph *9(0(2A—4y)-+ (r+ 1)(2r-+-3)a}% + 
2r+2 
+(2r+-3)wad+{2b w w (rapa) +241) 3a rap aro} 
2 arte 
/ + ba + 
+ ay" (r-+-1)h® Fs (r-+-2)h® egpenttig 
HABE) Het Der AF 
ar+2 


—(r+ 1)(2r-4-3)a*hAF2 4. (27-4 3) dF 
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+-gbA| hor¥(2r-+-2)08+ Fe p—ht# 8) _ 


(2r-+3)(r-+1)gah*-1(h| (y+ Fp F3] +09 =0 
(+) 


{—) 
Henceforth the following abbreviations will be used 
c = 4(3A°—4p"), d= 12(5A+6y), e = 24(4A?+ 9Au+ 6p?) 
f = 8(3342+-88Au +682), k= 9A+8p 
m = 2207+ 53Au+ 36p?, p = 39A?+-105Au + 80p? 
8 = d+40a, 


t = 3c— 
If r 0 


2e 
:. 3 


in (6.21) the following equations are obtained 
: ~~ 2 = 
abAs| htw—2h*w-+w| +A*{ch*w+adw|+0 = 0 


0 2 


6 0 4 2 
abA3! 2h®w— shtw+w! +A*{ehtw-+4-60a| bre —ah?| 


0 2 
—120abA| h*w— 30} 4+0=0 
. | eB e 
abA?} 3h w— 4hSw+w} +A?! fheu 4 


(6.24) 
6 2 
74 56a[2bw—3ah*|| a 


0 4 
—336abA| h*w— Su! +P = 0 
where 


(6.25) 
eee U 3 
pO = abd?! iinet a +4A[g[ kh? F?—3F?] 


‘° 2 6 
p® iit 2ns F3— 


pie de + 


2 4 5 1 
A{ mht ¥9—5| 20h F?— ph F*] + 9b] sh+ FP MFP g)) — 
Ls 
120ag{h?| %)p+F8jg+ F*] +-he) 
(+) 
0 


(6.27) 
ut = abA®| 3h°F*— 4h F834 F\ + 


2 6 7 1 
sA{ ph F3— 7|3a%h'F?—phF?] +-gb| 651% + Fp n° FP || “ 
336ag(h4] Pg + Fg 7] +e, 


(6.28) 
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2 
Eliminating Aw from (6.23), (6.24) yields the following equation: 
& £ % £, 
720abAw—3ash*A*w = abAsT + A*(th4w—120abw) + 


0 
+240abh27Aw+3h°Q—20, (6.29) 


0 4 6 


where T = —h®w+3h*w—2w. (6.30) 


Equation (6.29) may be = to ictaaio w from (6.23)-(6.25) save for 


a term in each involving Ate. The first two of these equations yield the 
same equation, viz. 


2 0 4 
LA%w = 2ab*A‘T + A*{b(3as-+2t)h4w+3ab(s—80b)w} + 
rS 
4-30 160ab2)h2A2w+- 2bA(3h2Q —2)+-380, (6.31) 
and (6.25) becomes 
2 ne 3 4 
RhtAtw = 4ab*h*A‘T +-A?{b(9as+- 4t)h8w+ Sabsw—480ab*h | + 
= 2 
+.A?{(3af+960ab*)h®w+ 336absw| + 
s & 
+-A/1008abs{ 5w—h*w) + 4bh4(3h7O—20)|+3s¥, (6.32) 


where L = 3a(48062—sd), R= 72a(7as+40b?). (6.33) 


Operating with A on (6.29) and eliminating A*w from the resulting equation 
by means of (6.31), the following equation is obtained relating the moments 
0 4 6 

w, w, and w: 


0 4 
(720ab—Sash*A)| 2ab*A4T+A*[b(3as-+ 2t)h4w-+ 3ab(s—80b)wo] a 
0 
+3(cs+ 160ab)h®A*w+ 2bA(3h70 —20)-+ 380 


0 4 0 ; 
== J abA‘T’ + A%(th*w—120abw) + 240abh*A*w+ A(3h*O—20)|. (6.34) 
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2 

Eliminating A’w between (6.31), (6.32) gives the following equation re- 
0 4 6 8 

lating w, w, w, and w: 


0 4 
Rh*{ 2ab*A4T +.A9|b(3as-+ 2t)h*w-+3ab(s—806)wo] + 


0 
+ 3(cs+ 160ab®)h*A%w + 2bA(3h70 — 20) + 380} 


0 8 4 
: L{ 4ab*h*AsT +-A*[6(9a8+ 4t)h®w+ 3absw —480ab*h%] +- 


4 
+ A2| (39f-+ 9¢0ab*)A* + 336absi] + Al 1008abs{5u—h'v) 4 
+ 4bh*(3h20 —20)| +38¥). (6.35) 


4 
The next step would be to eliminate w from (6.34), (6.35) but this is not 


necessary for the approximations considered below. 


7. Method of approximation 

It does not seem possible to give a method of approximation which will 
cover all cases simultaneously. It is for this reason that exact equations 
have been retained for as long as possible. The most suitable equation to 
use may frequently be determined by means of a preliminary solution 
using the crudest relevant approximation or may be obvious from known 
solutions to similar problems. It will be assumed below that units have 
been chosen so that |h| < 1. Equations (6.34), (6.35) will then give suffi- 
cient accuracy in the majority of cases. The equation 


0 
8u(A+ ph 3A*w+3(A+2y)p = 0 (p = —26%), (7.1) 
(-) 


for a plate under normal pressure p corresponds to that used in classical 


plate theory and will later be seen to be a first approximation in the 
0 


present theory. The functions h®A%w and p must depend in like manner 
on the parameter A. In the present case it will therefore be assumed 
that there exist constants M and w, independent of h, n, and v, such that 


Ane < Mh’-™ (n= 1,2,3,4; v = 0,4,6, 8), (7.2) 


where @ is the best possible value and is in general < 3. This condition 
arises from the assumption that Tien = O(h®). For suitably continuous func- 


(- 
tions, equation (7.2) would be satisfied if A"w (n = 1,2,3,4) are of like 
magnitude as far as h is concerned and 


\Aw| < Mh-*, (7.3) 
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If the latter equation is taken as the fundamental assumption it would 
have to be assumed that the function w is uniformly bounded for the 
ranges of values of x* concerned and is sufficiently smooth that differentia- 
tion with respect to z* does not change the order of magnitude with respect 
to h. Thus, functions such as sin(2*/h) and exp(a*/h) would have to be 
excluded because of the factor 1/h which appears on differentiation with 
respect to 2*. Functions of 2*/h which do not have this property would 
still be admissible. The displacement w may be a function of 2/h. From 


(7.1), and similarly from the equations obtained below, it can be seen that 
0 


0 

if the term A*w in the equation for w is to be retained and if the critical 

value of w in (7.2) is required when n = 2, then no significance can be 
0 

attached to an accuracy of the equation for w of order less than O(h*-”). 

Rejecting terms which in the light of (7.2) involve h7~”, (6.34) becomes 


0 4 0 
(240ab2—3es-+-td)h*A*w+ 120ab(6b—d)Aw+ 2400(3c-+a d)h2A2w-+ 

+ A{3(d—s)h?O — 2d®}-+- 72060 — 2bh?A?(3h7Q—2®) = 0. (7.4) 
When discussing the order of magnitude of the terms in an equation it 
will be assumed that the powers involved are h®, h',..., ie. the equation 
has been multiplied by a power of h sufficient to make the lowest power h°. 


It will also be assumed that f* + 0. If the latter function is zero it will 
(-) 


be found possible to divide the equations by h, or some power of h, and, 
in consequence, the accuracy of the resulting equations will be reduced by 
4 
this amount. The term involving A*w has to be eliminated from (7.4). 
When making the substitution from (6.35) terms which involve h*-” may 
4 
be omitted. To this accuracy A®w is given by the equation 


Rh*{2bA(3h°O — 20) +380} 


ne ee 
ws L{A[1008abs( Sw —h*w) + 4bh4(3h20 — 20) 436"), (7.5) 


or - H 
1680abLAw = 336ab Lh*Aw-+ 4abh*(84a+-d)A(3h2@ — 2) + Rh4O— LY. 

4 (7.6) 
Using (7.6) to eliminate A*w from (7.4) the following equation is obtained: 


0 
14Lh*{240ab*—3cs-+-td+24ab(6b—d)}A*w-+ 


0 
4 14.L{240b(3e-+ad)h®A*w—2A (60ah*O +d) +72060} + 
+-4ab(6b—d)(84a-+-d)h4A3(3h2@—2@)-+ 
+ A%(6b—d)( Rh*@— LY) — 28 Lbh?(3h?@ —20)} = 0. (7.7) 
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Equation (7.7) is the desired equation for w. It is accurate as far as terms 
of magnitude h°-”. To be precise, h’-” and higher powers have been 
rejected, but fractional indices are not considered here. The accuracy of 
(7.7) may be compared with the accuracy attainable with a fourth-order 
equation and that of the crudest possible theory as follows: 





Accuracy of (7.7) 
possible with a Accuracy of the 
Accuracy of (7.7) fourth-order equation crudest theory 





6—w 4-0 3-—aw 











The surface loading and body force terms ©, ®, ‘’ must be modified to 
the required degree of accuracy. This operation can be carried out without 
further assumptions in any particular case, but in order to proceed with 
the general theory it will be assumed that the body force and surface 
loading terms satisfy the following inequalities: 


(el <a a 


al, \F,) < Mh" (n= 0,1....), (7.9) 


and that these inequalities remain true after the operations A and A? on 
the functions on the left-hand side. It will be seen from (6.26)—(6.28) that 
under these conditions 
0 = O(h-), ® = O(h), ¥Y = O(A), (7.10) 
provided that #° + 0. If w = 3, the most accurate solution given above 
(—) 
is that with accuracy of order h’, viz. 
h2 S 
Su(A+ ayn A— i) 
3) 
; RS £ 
= 4(2A4 3u)h7A 1 —6(A+4 2y){h( 8g F8g4+ F%) + 3) + 
(-) X+) 


(3 


0 2 
+2(A+ 2p )APA 8 )94-(11A+ 12n)K9AF3—3MAF*, (7.11) 
(+) 


The corresponding fourth-order equation with accuracy of order h is 


0 0 
4u(d-+-p)h*Atw = 3(A+2y)/h( 289+ F) + . (7.12) 


(+ 


The crudest possible theory, viz. that with accuracy of order h®, is given 
by the equation 


0 
4u(A+p)h8A2w = 3(A+ 2p) #3. (7.13) 
(-) 
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Solutions to (7.12), (7.13) will in general satisfy the simpler basic 
assumption given in (7.3). In the case of (7.11), however, solutions of the 
forms A(h)exp(+2*v5/h) will arise and the value of w may exceed three, 
with a consequent lowering of the accuracy of (7.7). Even (7.2) is only 
a sufficient condition. The only real necessity is that the terms omitted 
from the exact equations (6.34), (6.35) are in fact negligible. As the value 
of w decreases, more and more of the body force and surface loading terms 
©, ®, ¥ become significant. 


8. Conclusion 


It has not been found possible to obtain a unique differential equation 
0 


for the mean value of the transverse displacement w which has a specified 


order of accuracy. The equation is seen to depend not only on the order 
0 


of accuracy required but also on the order of magnitude of w relative to 
the thickness of the plate. This, in turn, depends on the surface loading. 
However, for a wide class of problems, differential equations of specified 
orders of accuracy have been obtained. The accuracy attainable with a 
fourth-order equation and the significance of the surface loading and body 
force terms are clearly seen. The most important fact to arise from this 
analysis is that any attempt to refine the transverse displacement equation, 


e.g. by inclusion of a term of the form A #°, is likely to be accompanied by 
(—) 


a change in the order of the differential equation. 
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SUMMARY 


This paper discusses the effect of thermo-elastic coupling on the passage of the 
disturbance through an infinite medium with a spherical cavity when its tem- 
perature is suddenly changed. The hoop-stress at the cavity is seen to reach its 
steady-state value instantaneously; it then rises to a maximum, falls below the 
steady-state value, and approaches it again finally. The effect of coupling is seen 
to reduce this fluctuation of stress throughout and also to reduce the ‘shock’ 
character of the disturbance appreciably. 


1. Introduction 

THERMO-ELASTICITY is quite an old branch of study. However, steady- 
state solutions still form the major part of most investigations. Even in 
a number of cases of unsteady state a quasi-static treatment is made by 
rejecting the inertia terms, and solutions of the full dynamical problems 
are relatively few in number. 

In a majority of the earlier studies the equation governing the heat flow 
was taken to be that given by Fourier. This classical equation assumes 
that the heat conduction is independent of the elastic state of the body. 
This hypothesis was recently re-examined and the equation reformulated 
(1, 2) by taking account of the effect of rate of strain on the flow of heat. 
In the classical case the equation of heat flow can be solved independently 
and the equations of elasticity solved subsequently. In the modified 
equations this is not possible; the equations of elasticity and heat con- 
duction are coupled and have to be solved simultaneously. This increases 
the complexity of the problem; so it would be interesting to study how 
far the modifications affect the conventional solution. 

The present discussion concerns such an investigation. In particular 
we consider an infinite medium with a spherical cavity. The surface of 
the cavity is maintained stress-free, but its temperature is instantaneously 
changed to some different constant value. The ensuing disturbance is 
the object of our study. The full dynamical problem, with the classical 
heat conduction equation, was recently discussed by Sternberg and 
Chakravorty (3). We refer to this as the conventional solution. 

(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961) 
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The present study is, of course, valid only within the framework of the 
classical linear theory of elasticity. This limits the range of allowable 
changes in the temperature. Further, we neglect, as a simplification, the 
variations with the temperature in all the material constants. The basic 
equations are those derived by Biot (1). 

The coupled equations are solved by the use of the Laplace transform 
technique, the solution of the general problem being obtained in the form 
of inverse transforms. The hoop-stress is studied in detail. In such 
dynamical problems the problem of inversion is often very involved since 
it is not always easy to obtain expressious which are amenable to calcula- 
tion. This is so for this problem, even in the conventional case—and the 
coupled problem presents greater difficulties. Thus a series development 
valid for small values of the time is given first. However, this series is 
of little use to get an overall picture of the fluctuations of the stress, 
and it is this which would help us to understand the nature of the effect 
of coupling. Since the thermo-elastic coupling constant is found to be 
small, a perturbation solution for the hoop-stress at the wall of the cavity 
can be obtained. This enables us to obtain an integral representation for 
the solution in terms of tabulated functions only. We feel that, whenever 
inversion is possible in the conventional case, this method of perturbation 
enables us to obtain the solution for the modified case as well, the precision 
achieved being the same, to the order of the perturbation developed. We 
have used this perturbed solution to obtain a power series solution valid 
for those ranges of time which are of major interest. 

The instantaneous change in the temperature of the wall of the cavity 
gives rise to thermal and elastic disturbances which travel through the 
medium. With regard to the circumferential stress, it is seen that the 
disturbance consists of two parts as in the conventional case (3). The first 
part is diffusive in character and the second part travels as a wave. At 
the head of the front there is a discontinuity in the stress, which explains 
the use of the term ‘shock’. The effect of coupling is seen to reduce the 
strength of this shock to a very large extent, except at points very close 
to the cavity. As in the conventional case (3), the hoop-stress at the 
cavity attains its steady-state value instantaneously, rises to a maximum, 
and then falls below it, approaching it finally from below. The effect of 
coupling is seen to reduce this deviation from the steady-state value 
throughout. 


2. Basic equations 
We consider an infinite isotropic medium with a spherical cavity of 
radius a. Let (7,, 7), Ty) be the radial and the circumferential stresses. 
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Assuming complete spherical symmetry the radial displacement u is the 
only displacement. We then have 


ou u 
T = Ty, oy Oe ae om OS (1) 


where ¢,, ég, and e4 are the radial and the circumferential strains. 
The stress-strain relations are 
T, = Ae+ Qyve,—a(3A+ 2y)T, 
T) = T; = Ae+- 2peg—a(3A+ 2)7, (2) 
where e = e,+-e9+e,, A, » are Lamé constants, « is the coefficient of linear 
expansion, and 7’ is the excess of temperature over the absolute tempera- 


ture 7, when the medium is unstressed. 
The equation of heat conduction now takes the form 


ry2 = or ee 
KVT = 6, + Tya( B+ 2p), (3) 


where V2 = @/dr?+-(2/r)(@/ er). 
Here K is the thermal conductivity, c, is the thermal capacity at constant 
volume, and ¢ is the time. 
Finally, we have the equation of equilibrium 
oT, 2 ou 
—'+4+—(T,—T}) = pp—» 4 
or + ms, r a) Po ot2 ( ) 
where py, is the density. 
We reduce the above equations to non-dimensional form by the follow- 
ing substitutions: 


1 
TR; 9%, 04) = = (T,, T), Ty) 
(op; %, % on 1% 


(u,r) = a(U, RP), T = 7,9, 
oe A 
Po 
Then the equations (2), (3), (4) assume the forms 
_ OU , A[eU , VU) $a (3A+2p), 
or = oR oleat a ea 
U  Afau , 2U\  aT,(3A+- 2p) 
phat! ae aieneenmnis 
+ sen 4 2 . 
C,, av 9 aav(3A-+-2) de 
K @r K @r’ 
fen 3 Pov Hu 


oR ROR) = on oe 


{= —T; v? 
v 


og > @ 


eee 
V0 = 
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Using (2a) in (4a), we finally obtain 
DD,U—a, DO = 


fa 0) 
D, DO —24— = ag 


aoe ( BA+ 2) 
. A+ 2p 


are non-dimensional constants and 


> 


ao © 
om oy 
3. Solution of the problem 
Taking Laplace transforms of (5) and noting that U = 0, 0 = 0 for 
t= 0, we get (DD,—s*)U = 4 DO \ 
(D,D—a8)0 = a58D,0)’ 
where the bars denote the transforms of the corresponding functions 
throughout. 


Operating by DD, and D,D respectively on these equations and using 
the remaining one, we finally obtain 


{L?—(m3-+m3)L+m}m3}0 = 0) 

{M?—(m3+-m3)M +m? m3}6 = 0)’ 
where L = DD,, M = D,D are two operators, and mj? and m3 are the 
roots of the quadratic in m?, given by 

m*—{s*-+-(1-+-€)a, 8}+ a, 8* == @, 


(6) 


(7) 


with € = a, a;/a,, denoting the thermo-elastic coupling constant. 

We thus note that the temperature and displacement depend on two 
non-dimensional constants, « and a,, the former depending only on the 
material constants and the latter on the linear dimension as well. 

We frequently refer to the orders of magnitudes of these quantities 
though we do not use them very much. As an example we give below 
their values for aluminium, taken from (4). 


v = 6-30x10° cm/s; TT, = 300° K, 
K/c, = 0-4552 cm*/s;  v = }, 
my = 13-84x108a = 107, =e = 00-0269. 


The two fourth-order equations (7) possess, in general, four independent 
solutions each, but since U and © must vanish at infinity, we exclude 
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two solutions of each. The solutions of the first of (7) are modified 
Bessel functions of order 3 and can be expressed in terms of exponentials. 
Thus we have 


(8) 


where the a; and 6; are functions of s. Since the expressions (8) must 
satisfy (6), only two of these functions are arbitrary. Also since the 
solutions must vanish at infinity, the square roots are always so chosen 
that they have positive real parts. 

Now the boundary conditions of our problem are 


p= 90 (for R = 1, for 7 > 0), 
T =T,= xT, (for R = 1 and for r > 0), 


where y is a constant. These give 


" (for R = 1). 


oR 


b, os: 6 
i » 
Oy mM; 


a,e~™ ‘a —a,e-™ em 
my(%48*+-2+-2m,) — mo(ags?+-2+2m,) = 8A’ 








where ig ae: bing 


A = (m,—mzy)[ (x, 8?+ 2)(m,+m,)+ 2(s*+-m, m,)]. 


Using these values we can obtain all the relevant quantities. To see the 
effect of coupling, we study the circumferential stress in greater detail. 
Its value is 


lip sn pk [ (x, 8?R?@—1—m, R—m? R*)(a,8?-+ 2+ 2m,)exp(—m, R,)— 


— (a8? R?—1—m, R—m§ R*)(a,8*+ 2+ 2m, )exp(—m, R,)], (9) 
where R, = R—1. 
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4. Calculation for small values of the time 


The above expression, as it stands, is too complicated for inversion. 
However, we can obtain an estimate, valid for small time, by expanding 
it in inverse powers of s. 

This gives 

Boa. Sit exp{—R,y(«28)} x 

Oy X R 





x [as ech Ahead ag Ri+(1—e)ag a4—a4}. soe 


_ exp{— meal) 
R 8 


wade 
+5{(4—)((0—o ; e(4 v2" l 


On inversion, we obtain 


1 Rw Rv. 
= plssert a+ a, €x§ R, v7 erfe a+ 


+-ag{fe*a§ a, RE+(1—e)ay—1} x 


x (ere Ri erfe =a Ry |(Jexp(— “icy).. ea 





a 


pexP(— jea, R,)H(7— Ry) (a 1)—(r—R,) x 


x {(a—0(i4—oeat R,+=—(1—)a) +pt cag}. (11) 


The solution corresponds to expansion up to terms of order s~? only. 
The first part clearly represents the diffusive part and the second represents 
a stress wave travelling with a velocity corresponding to that of the 
dilatational wave. This is analogous to the conventional case (3). For the 
stress-discontinuity across the wave front, we obtain 


%9 


(~24) -(-24) a 8! ox —Je0, R,) 
% X/r=R:-0 % X/r=Ri+0 R 


ASS. 4 
= [om Ror — te R,). (12) 


It is this discontinuity of the stress at the head of a wave that gives it 
a ‘shock-like’ character. In the conventional case (3) the exponential 
term is absent and the strength of the shock then varies inversely as the 
radial distance. The decay here, due to the exponential factor, is far more 
rapid. From the numerical values quoted, we see that the quantity «a, 





ON SPHERICALLY SYMMETRIC THERMAL SHOCK 81 


is of the order of 105, So, except at points very close to the cavity, the 
jump is negligible. The expansion (10), from which (12) is obtained, is 
of course valid for small time only. However, we note that this jump is 
due solely to the first term of the second part of the expression (11). 
Higher powers of s~! would contribute nothing to this jump. So, though 
not strictly proved, it is plausible to infer that this law of decay persists 
for a longer time. 

To obtain the steady-state value, we note that the singularity of (9), 
which has the largest real part, is s = 0. Expanding (9) in ascending 
powers of s, we obtain 

ome = se gat p) + higher powers. 
Thus aS re sla 2) as 7 -> 00. (13) 


mx 2 
From (13) we see that the steady-state value is unaffected by the thermo- 
elastic coupling constant, as expected. 


5. Calculations for R = 1 
In this case we have 


Gp if (3a4—2)s? } 





8 (xg 8?-+-2)(m,+m,)+ 2(s?-+-m, m,) 


On expanding in ascending powers of s, we get 
a7) 1 (3a,—2) 


i 


—— = —- +4 —_*__—.. g$-+-higher powers, 
a4 X 8 2,/{a9(1+-€)} 8 po 


and inversion then gives, for large values of the time 
Rs a Mee pap ES ae 
%4 X 21(3)ytae(1+-«)} 
From this we see that the magnitude of the steady-state value is attained 
from below, as in the conventional case (3). We subtract the steady-state 
value and let the difference of og and its steady-state value be of. We 
further introduce, for convenience, a quantity o proportional to this o}, 





ay X(Bag—2) (xg. 8? ++- 2)(mM +- My) + 2(8?-+- Mm, My) 
Expanding this in inverse powers of s, we obtain 
-___ [laf , (1—he)ag (1—Be)ah  (1—3e+Fe*)ag j 
= “| ait Bens AO st x 


o= 





(15) 





(16) 


where, because of its order of magnitude, we have neglected terms that 
are divided by ay. 


5002 .53 
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Using the numerical value of «, quoted above we find 
o = 7,—0°75257} + 0-493373—0-22947] +-0-153373—..., 
where 7, = a 97. 
The following table gives values derived from the above. 
Th ° or o2 o3 o4 ors 
oe ° 01026 01829 02136 o2752 03359 
We see from the above power-series development that the effect of «, 
the thermo-elastic coupling constant, is to reduce the stress-deviation, 
as represented by oc. However, since a, is of order 10’, this gives values 
for very small r, of the order of 10-’. In order to get an estimate for larger 
values of r, we develop another expression for 6, obtained by perturbation 
of «, a small quantity. 
Expanding ¢ in powers of «, up to its first power, we obtain 








eS. Vs __ El v(s?-+- 2p) 
ag  (V8-+-Varg)(8+-8,)(8+8,) 2 (V8-+Vavg)®(8-+8,)*(8+-82)? 
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wae bi2—=ptig, p=Nay 9 = V(2q—1)/x. 
After a partial fraction decomposition, we obtain 
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We note that whenever the inversion of F,(s) can be achieved in terms 
of tabulated functions, so can that of F,(s) by using some modifications. 
The expression F,(s) corresponds to the conventional case. In the present 
case F,(s) cannot be inverted to give a result in terms of tabulated func- 
tions. Using the decomposition 


1 1 Note Vote , 8» 


a 


8+8, Va(ve+Vay) | V8(8+8,) 








(W8+WVag)(8+8)  a+8, 
we can express F,(r) in terms of error functions \ ith the complex argument 
p+iq. However, since these are not tabulated, the solution does not offer 
much advantage in performing calculations. So we express these in 
alternative but equivalent forms by the use of the convolution theorem. 
Though these are in an integral form, we can get an approximate value 
of o for ranges of r which are of major interest. 
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By the use of the known formulae, we get 
L-'F,(s) 


1 . Qe J e ' 1 
a ao? —P(T-%) gin —u)| e™"“ erfe./(a. u - > du 
= sin gr + Ps i sin q(r | (aU) Ina, %) 


Tt 


= | e™“ erfc Maguye**-| 008 q(r—u) —{Gsin ar—w)| du (18) 


after integrating by parts. Thus 
: rd ro v p. fv 
xP (r) = | e’erfe vv e”"'™-! cos q|——r7)+~sing|——r}|dv, (18a) 
3 Xe q Xe 
0 
where we have taken v = a, uv. Similarly, 
L-"F,(8) 
Pp f 
g i 
0 


K | (1+ Sag u-t 2nd wt)ere" erfe V(ag u)—2(2+ Xo 0 /(=)| du. (19) 


[2sin g(r—u)—gq{p cos q(r—u)-+-q sin g(r —u)}(r—u) Je-"-” x 


Again, after integrating by parts, we obtain 


2F(r) = | (20+ we ef sv—20 /(2)] x 
7 
0 


moe | 2 sin g(u—r)+q(u—r) qsin q(u—r)—pcosq(u—r)}}em=-") dv. 
du ¢ } 


, (19a) 

In order to evaluate these expressions, we expand the trigonometric 

and exponential functions in power series. The range of validity of this 

series is the same as that of the exponential series and the latter is known 

to be convergent for all values of the argument. However, we do not use 

this series beyond + = 2-5. We then integrate (18a) and (19a) term by 
term. For this we require the following integrals: 


[ (2—*)"ererfe vv dv = 2 /(*2") _ el Mal 
J \ag 1 )3.5...(2n+-1) 
0 


Xer v n Dy2 vr los. Dae v — AF _ — 2r)"n! 
| (— 7} (ze +-v)e" erfe vv 2» /(2)) dv = 2 /| 2 \y pea (20) 
0 


The values of these integrals are approximate. Terms of order unity or 
lower are neglected as compared to va,. We further use the numerical 
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Using the numerical value of «, quoted above we find 
o = 7,—0-7525r} + 0-493372—0-22947}-+-0-153373—..., 
where 7, = a rT. 
The following table gives values derived from the above. 
% ° or o2 o3 o4 os 
o ° 01026 01829 02136 02752 03359 
We see from the above power-series development that the effect of «, 
the thermo-elastic coupling constant, is to reduce the stress-deviation, 
as represented by o. However, since a, is of order 10’, this gives values 
for very small 7, of the order of 10-’. In order to get an estimate for larger 
values of r, we develop another expression for ¢, obtained by perturbation 
of «, a small quantity. 
Expanding ¢ in powers of ¢, up to its first power, we obtain 
ee. vs Re vs(s?+- 2p) 
ay  (V8-+Varg)(8+8,)(8+8,) 2 (V8-+Vag)*(8-+8,)*(8-+8)" 
= F,(s)—}ea, F,(s), say, (17) 
where S2=ptiq, p= liam q = V(2ag—1)/ ax. 
After a partial fraction decomposition, we obtain 
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We note that whenever the inversion of F,(s) can be achieved in terms 
of tabulated functions, so can that of F,(s) by using some modifications. 
The expression F,(s) corresponds to the conventional case. In the present 
case F,(s) cannot be inverted to give a result in terms of tabulated func- 
tions. Using the decomposition 


1 SMe ak f Vag 8 

(V8+~Va)(8+8,) al Va(V8+- Vay) | aa | 
we can express F,(r) in terms of error functions with the complex argument 
p-+iq. However, since these are not tabulated, the solution does not offer 
much advantage in performing calculations. So we express these in 
alternative but equivalent forms by the use of the convolution theorem. 
Though these are in an integral form, we can get an approximate value’ 
of o for ranges of + which are of major interest. 
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By the use of the known formulae, we get 
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= | e*“erfc,/(a, uje-#r-»| 008 q(t —u) —{sin atr—w)| du (18) 
after integrating by parts. Thus 


12F,(r) = | ererte vvenae-r| cong( 1) +2sin (=—-)| dv, (18a) 
0 


Xe Xe 
where we have taken v = a,u. Similarly, 
L-"'F,(8) 
T 
pf 
q° 


0 


m | (1-+ 5g u-+ 208 wthemenfe (ag u)—2(2+ aw), [(***) du. (19) 


[2 sin g(r—u)—gq{p cos g(r —u)+-q sin g(7—u)}(r—u) Je-?7-” x 


Again, after integrating by parts, we obtain 


XeT 


= [ (208+) erfe vv—20 /(2)] x 
é 7 
0 


& 7, (2sin q(u—r)+-q(u—r){qsin g(u—r)—p cos q(u—z)}}e*-”) dv. 


* \du @ 
(19a) 
In order to evaluate these expressions, we expand the trigonometric 
and exponential functions in power series. The range of validity of this 
series is the same as that of the exponential series and the latter is known 
to be convergent for ail values of the argument. However, we do not use 
this series beyond + = 2-5. We then integrate (18a) and (19a) term by 
term. For this we require the following integrals: 


,, n _. De\Ny t 
[(2- e’erfe vv dv = 2 | ms, os... 
J} \as a )3.5...(2n+1) 
0 


XeT a ne a sey v a4 %eT\ (—2r)"n! r 
| 5 *} [ize + ve" erfe vv 2» /(2)] dv = 2/227) ti)’ (20) 
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The values of these integrals are approximate. Terms of order unity or 
lower are neglected as compared to va,. We further use the numerical 
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values p = 4 = cos(}7), q = sin(47), corresponding to vy = 4. We then 
obtain 


7, _» < fsin}(n+1)r a.) (—2r)" 
Jz oe Mier “5 | Ee ~ Pee 


= fi(r)—4ef.(7), say, (21) 


where a, = 3 gl2sin }(n+1)m—-(n-+1)sin }ar cos }(n+1)m]. 


Since ¢ is very small, the valuc of o is very little modified. In order to see 
the relative magnitudes of the two functions f,(r) and f,(r), we tabulate 
them separately below: 


T | ° ool or o2 o-4 os 06 





Sit) ° o-2108 0° 3880 0°4660 04780 04755 
S,{r) ° 02757 03375 0° 3483 03278 03026 


T I's 2 25 





Sir) o°2140 0°0440 —oogIg 
}a(r) 0°0397 —o°oosi — 

From the table we see that f,(7) passes through zero at + = 2-2 approxi- 

mately. Though the modification caused by « is quite small, its effect is 

to delay the time when o passes through the value zero. Also, in the 

range considered, it tends to reduce o. 


6. Conclusion 

The above analysis leads us to two important inferences about the effect 
of thermo-elastic coupling: (1) the stress-discontinuity is reduced to an 
almost negligible amount except at points close to the cavity; (ii) fluctua- 
tion of the hoop-stress at the cavity about its steady-state value is reduced 
throughout. This, of course, assumes that + = 2 is large enough in the 
history of the variation of the stress. Also the perturbation method 
developed in the paper simplifies considerably the problem of inversion 
of transforms. 
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SUMMARY 

A general method is presented for the solution of the differential equations for 
the stability of a deflected and rotated bar subjected to variable axial loading and 
constrained elastically along its length against deflexion and rotation. 

The method is used for the study of the stability of compression chords for 
through-bridges. By using Fourier sine series for the deflexion and rotation curves 
of the chords general solutions are obtained for the most important loading condi- 
tions encountered in practice. 

The convergence of the series solutions and the practical application of the solu- 
tions for the design of compression chords for through-bridges are discussed. It is 
concluded that the analyses presented in the paper should, within the limits out- 
lined, be applicable to the design of most of the usual types of through-bridges 
without too elaborate computation. 


1. Introduction 

THROUGH-BRIDGES are bridges where the cross-girders are placed below the 
level of the compression chords of the main girders, and where there is 
no lateral bracing between these chords. In bridges of this type, the cross- 
girders and the vertical stiffeners of the main girders will form U-frames, 
which will support the compression chords in such a way as to resist any 
deflexions and rotations of these chords. The ability of the chords to 
withstand compressive loads will, therefore, not only depend on the 
flexural, torsional, and warping rigidities of the chords themselves, but 
also to a great extent on the stiffness of the U-frames. 


i.1. A great amount of work has been carried out in the study of the 
behaviour of the compression chords for through-bridges, and an account 
of the research in this field during the last seventy-five years has been 
given by F. Bleich (1). Most of the design methods proposed are, however, 
very cumbersome for practical use and most of them are based on the 
assumption that the chords remain straight up to the point of buckling. 
In practice, the live load will usually be applied on to the flexible cross- 
girders which will deflect and cause the tops of the U-frames to deflect 
and rotate. In most types of through-bridges the ends of the compression 
chords will be prevented, to a certain extent, from deflecting laterally 
and from rotating in the planes of the end U-frames. The application of 

(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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live load on to the cross-girders will, therefore, lead to deflexions, rota- 
tions, and stresses being imposed on the chords which clearly will affect 
their ability to carry compressive loading. 


1.2. In this paper a general method is presented for dealing with the 
behaviour of the compression chords for through-bridges. It is based on 
the energy method which was first used by Timoshenko for the study of 
problems of this type (2). The U-frames are assumed to be replaced by 
a uniform elastic medium, and the axial loading in the chords is expressed 
as a continuous function of the lengths of the chords. The effect on the 
stability of the chords produced by placing the live load on to the flexible 
cross-girders is considered. This implies that the chords will deflect and 
rotate long before the buckling loads for the straight chords are reached. 


2. The general method 

Consider a straight bar subjected to axial loading and supported laterally 
along its whole length by a uniform elastic medium. At a certain value 
of the axial loading, the bar will change suddenly from the straight form 
to a deflected configuration without any change in the value of the axial 
loading. Using the compressed but undeflected state of the bar as reference 
position, the following energy equation can be set up: 


V.tV,—W =0, (1) 


where during the buckling of the bar, V, is the work done in the bending 
and the twisting, V, is the work done on the elastic medium, and W is 
the work done by the axial loading. 


2.1. The expressions for V,, V,, and W will all be integrals, extending 
over the whole length of the bar and containing functions of the deflected 
and rotated form of the bar. W will also contain the expression for the 
axial loading, which, normally in practice, can be expressed as a continuous 
function of the length of the bar. 


2.2. Expressing the actual deflexion and rotation curves of the buckled 
bar by series which satisfy the boundary conditions, and introducing these 
series into equation (1), this becomes after integration: 


F,+F,—PF, = 0, (2) 


where F,, F,, and F, contain functions of the parameters a, for the de- 
flexion curve and b, for the rotation curve of the buckled bar. P is the 
value of the axial loading at a chosen section of the bar. From equation 
(2) the value of P can be expressed as 


P= (4, +,)/F. (3) 
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2.3. The coefficients a, and 6, are parameters which must be determined 
so as to give a minimum value for P. This condition can be expressed in 


the following form: 

éP 1/eF, ,@F, F4+F, | 

—— = ei! + — 3 1 3 — 3] — 0 

Ca, F; oa, éa, F; oa, 

aP _ leh, a K+h eh) _ 

ob,  F;,\eb,, ° ab, yy &] . 
Since F, ~ 0 and (F,+F,)/F, = P from (3), these equations become 
ph 


éa,, 


oF, eF, 
a a os 0 
scooelt <p 


n Y n ‘ (5) 
oF, , oF, eF, 
0b, ob, ob, 
By giving n in the above equations the values 1, 2, 3,..., 7,..., two sets of 
an infinite number of linear homogeneous equations are obtained. For 
buckling, the coefficients a, and 6, must be different from zero, and the 
buckling value for P can consequently be obtained by equating the 


determinant of the coefficients a,, or b,, to zero. 


2.4. The buckling value for P can also be obtained by solving the 
differential equation for the buckled bar. This differential equation can 
be obtained by either of the following two methods: 

(a) by considering the equilibrium of an element of the buckled bar, or 

(6) from the energy equation (1). 

Method (6) above is an extremum problem which can be reduced to 
the solution of the Eulerian differential equations with regard to the 
particular boundary conditions for the case under consideration. 


2.5. The energy equation (1) consists of integrals which are of the form 
L 
= | F(z,u,¢,u',d',u", 6") de, (6) 


where primes denote differentiation with respect to z, an abscissa measured 
along the length of the bar. It is required to find functions u = u(z) and 
¢ = (2) which will make the expression (6) for J stationary. This can be 
achieved by the solution of the following Eulerian differential equations: 


éF d(éeF\  d*(éF ' 
— — —{ — —-| —— } = ( 
ou aa ‘ Ane) . 


eF =Asg sls) = ° 


(7) 
a dalag’) * daag 


These differential equations are identical with those which are obtained 
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by using method (a) mentioned above. Provided the exact forms of the 
deflexion and rotation curves are used, results obtained by the solution 
of equations (5) or (7) will consequently be identical. 


2.6. In the case when the bar does not remain straight up to the 
point of buckling, buckling will not take place suddenly and the energy equa- 
tion (1) cannot be set up so readily. Instead a variational equation will 
be considered. This variational equation will be similar to equation (2), 
and is chosen so that, when substituted in the Eulerian equations (7), 
differential equations identical to those obtained by considering the equi- 
librium of an element of the deflected and rotated bar, are obtained. 
Representing the actual deflexion and rotation curves of the bar by series 
which satisfy the boundary conditions, and performing the integration, 
the following general expression for equation (2) is obtained: 


F,+.. +F,— PF, i+ rte +f) = “a 0, (8) 


where F,,..., F,,..., F, contain functions of the parameters a, for the 
deflexion curve and 6, for the rotation curve of the bar. These coefficients 
must be determined in such a manner as to give a minimum value for P. 


By comparing with equations (4) and (5), this condition leads to the 
following two equations: 


6a, n 


oF, + Oh p(n, 49H) _ of} 
&,  & By 

By giving n in the above equations the values 1, 2, 3,..., n,..., two sets of 

linear non-homogeneous equations are obtained. The deflexion and rota- 


tion curves can, for any particular value of P, be determined by the 
solution of these equations. 


ia, 


ea, aa, 


(9) 


n 


2.7. The method outlined above for the solution of stability problems 
is of value in the study of the behaviour of the compression chords for 
through-bridges, since, by using Fourier sine series for the deflexion and 
the rotation curves, very accurate results are obtained by using only a 
few terms of the series. Live load covering the whole span of a bridge 
will cause deflexion and rotation curves which are always symmetrical 
about the midspan, and thus only odd values for the Fourier coefficients 
are used. For the even values of the coefficients, the linear equations will 
become homogeneous and thus correspond to the buckling of the straight 
chords. For the particular case when the cross-girders are infinitely stiff, 
the chords will remain straight up to the point of buckling, and the linear 
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equations will be homogeneous both for even and odd values of the 
coefficients. 


3. Defiexion and rotation of the compression chords due to sym- 

metrical loading on the cross-girders 

In this paper only through-bridges having two parallel main girders of 
uniform depth are considered. The cross-girders are assumed to be evenly 
spaced along the whole span, placed at right angles to the main girders 
and rigidly connected to the vertical stiffeners. The bridge is assumed to 
have such rigidity in the horizontal plane as will prevent the main girders, 
at any section, from moving sideways in the same direction. The top 
chords are assumed to be straight before any live load is applied on to the 
bridge. Only a uniformly distributed strip load covering the whole span 
is considered, and the resulting bending moments and axial loads in the 
chords are assumed to vary sinusoidally. Only the elastic behaviour of 
the chords and the U-frames is considered. The flexural and torsional 
rigidities of the tops of the U-frames, in the plane of the frames, are 
assumed to be replaced by an equivalent elastic medium. The torsional 
rigidity of the vertical stiffeners, in the horizontal plane, is neglected. 


3.1. Consider a unit length of the elastic medium which supports the 


compression chords in the lateral direction. The medium is of the same 
form as the actual U-frames, and the horizontal and in what follows the 
vertical members are referred to as the cross-girders and the vertical 
stiffeners respectively. 


3.2. Due to the symmetrical loading on the cross-girders, these will 
deflect and cause the tops of the vertical stiffeners to deflect and rotate 
through A and ® which, by reference to Fig. 1 are given by 

A= 6 | 


(10 
o—0@) (19) 


Fie. | 


Unless, then, the actual deflexion and rotation of a chord at a section are 
A and © respectively, it will be acted upon by reactions from the elastic 


repesnee 


Rear er cove Sees 
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medium. When the actual deflexion is u and rotation ¢, as shown in Fig. 2, 
these reactions R and M will be proportional to @ = A—u and ¢ = ©—¢ 


respectively. Using the principle of superposition, i and ¢ can, with 
reference to Fig. 2, be expressed as follows: 


a = Rk+Mlky | 
A 


ll 
$ = Riko +M/ky (11) 





4 
uu 


“ 








4 
¢ 











where the properties of the elastic medium are defined by: 


k = equal forces acting at A and B to cause unit deflexions at A 
and B; 
k, = equal moments acting at A and B to cause unit rotations at A 
and B; 
ko, = equal forces/moments acting at A and B to cause unit rotations, 
deflexions at A and B. 


From (11) the réactions can be expressed in the form 
R = £,(A—u)—£,(®—¢) 
M = £(®—4¢)—£,(A—u) J’ 
E, = kkhy/(kii— kok), 
§, = kkgkgy/(kii—kyk), 
E, = ky kg, /(kii— Kok). 
3.3. The differential equations for the deflexion and the rotation of a 


chord \vhen subjected to an axial loading P, = Psinzz/L = Psinwz, can 
be written in the form (3) 


where 


Elu'¥ + Pu’ sin w+ Pwu' cos wz— R+ Py, ¢" sin wz+ Py, wd’ cos wz = ml 
EV ¢'*—GK¢"+ Pri’ sinwz+ Pr? wd’ cos wz—M + Py,u” sin wz+ 

+ Py, wu’ cos wz = i. 

(13) 
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where, for the chord,. HI is the flexural rigidity, GX is the torsional 

rigidity, ET is the warping rigidity, r, is the radius of gyration with 

ect to the shear centre and y, is the distance from the shear centre to 
ntroid. 


The solution of equations (13) is obtained by the method outlined 
beginning of the paper. The following variational equation is 
ed: 


L 
- ET $"*+GK$)de+ | (E,u?-+£5¢*—26, ud) d2— 
0 


L L 
-P | sin w2(u'?-+-4'*+ 2ygu'p’)dz—2(£, A—£,®) | u dz+ 
0 D) 

I 


+2(f,A—£,0) [ ¢dz = 0. (14) 


0 
vtional equation wil!, when substituted in the Eulerian differential 
is (7), give the differential equations (13). Assuming the deflexion 
tion curves to be represented by the series 


u = >a, sinwnz | 
¢ = > b, sinwnz \ 


e required solution of (13) will be identical with the solution of (14) 
ed in accordance with (9). Proceeding in this way, the solution 
tions (13) can be expressed as follows: 


EL Iw®rn* +4, L—2wPB,|+2wP > a,B,.— 
bal dE, L+2wPyoB,]+2uP yo > beBns = 2(E, A—Ey)/wn, 
8 
bf, L+ 2w PyyB,|+ 2w Py, > 4,By.+ 
8 
-b, [4 BT want +46 Kwan? + b€,L—2w Pri 8,,|4+2wPré > 6,8, 
8 
= 2(€,2—€, A)/wn 
(oe == 1, 3,6....;.¢ = 1,3, 5....,; @ #8) (16) 
B,, = n*(2n?—1)/(4n?—1) 
Bas = ns(n*+-s*—1)/(n*+-s*— 2n?s?— 2n?— 2s*+-1), 
3.5. For the case when the chords have negligible torsional and warping 
rigidities, the equations (16) simplify to 
a,| 4B Iw*rnt +4kL—2wPB,|+2w0P ¥ a, 8, = 2kA/wn 
8 
(xs == 1,3, 5,...; @ = 1,3, 5,...; 2 #8). (17) 


(15) 
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4. Defiexion and rotation of the compression chords due to un- 
symmetrical loading on the cross-girders and lateral wind load - 
ing on the main girders 
In order to determine for which deflexions and rotations of the chords 

the reactions from the elastic medium vanish, the following three different 

conditions will be treated separately, considering a unit length of the 
elastic medium. 


4.1. Effect of the unsymmetrical loading 
Due to the unsymmetrical loading on the cross-girders, the deflexions 
and rotations at the tops of the U-frame will, with reference to Fig. 3, 
be given by: bse: Agi oe | 


(18) 
O4,= 04; 9p, = On 


Ap, Sar 
<«_—— 


Fic. 3 

4.2 Effect of the difference in the main girder deflexions 

Due to the unsymmetrical loading on the cross-girders, the main girders 
will be subjected to different loading intensities g, and gp. The difference 
between the main girder deflexions at a section z can be expressed as 

2 = N4— Np = (Ga—Va)(2*— 2208 + L*2)/24 EI, (19) 
where EI is the flexural rigidity of the main girders. The resulting 
deflexions and rotations at the tops of the U-frame can, with reference 
to Fig. 4, be expressed as follows: 
A, = en/d = 2Q(z*—2 128+ Lz) 


ge F (20) 
®, = n/d = Q(e*—2L2°+ L*2)/e 
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4.3. Effect of wind loading on the main girders 
With wind loading as shown in Fig. 5, the resulting deflexions and 
rotations of the tops of the U-frames can be expressed as 


Maw = (Wy +Awg)/ky 
Anw = (wgt+Auy)/k, 
Diy = (Wy t+yep)/hs 
Dpy = (wet+yus)/ks 


Asw Aaw 
B 


/ 











where A = k,/k, and y = k,/k,, and the properties of the elastic medium 
are defined by: k, = force acting at A to cause unit deflexion at A; 
k, = force acting at A to cause unit deflexion at B; k, = force/moment 
acting at A to cause unit rotation/deflexion at A; k, = force/moment 
acting at A to cause unit rotation/deflexion at B. 


4.4. Using the values obtained above, the deflexions and rotations of 
chords for zero reactions from the elastic medium are: 


Ag = 44i—4,+Aawi Ap = Apet+A,+Aaw | 


: (22) 
D, = O,,—9, 49047; Of = 9p,+9,+9py 


When at a section the actual deflexions are u, and u, and the rotations 
are ¢, and ¢p, the corresponding deflexions and rotations from the posi- 
tions of zero reactions from the elastic medium are, with reference to 
Fig. 6, given by 

ty = Ay—uy; ty = Ap—uz | 


& = D,—¢; op = Oz— bp \ 


The resulting reactions from the elastic medium can, using the principle 
of superposition, be expressed by the following equations: 


(23) 


k, Ry +-Aky Rp+k, My+yk, Mp = ky ks ty 
Ak, Ry +k, Rp+yk, My +k, Mp = ky at 
ky Ryt+-yky Ry +key M,+ pkg My = ky kd, | 
yk, Ry t+k, Ry+pky My+ky My = kghy dy 
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where p = k,/k, and k, = moment acting at A to cause unit rotation at 
A and k, = moment acting at A to cause unit rotation at B. 


4, 


[wp ty! 
ee 


Rz Ms 


B i\ / 


%, 











4.5. The following assumptions are now made: 
Up = Py = BU; op = Udy = po- (25) 
The correctness of these assumptions will be discussed later. Using (23) 
and (25) the reactions due to the elastic medium can be expressed by 
means of the solution of (24), as: 


R, = —(A—pB)d—(C—pD)u+-A®,— BO,+CA,—DAz 
Ry = (B—pA)$+(D—pC)u— BO,+ A9,—DA,+CAz 
M, = —(F—pH)¢—(A—pB)u+ FO,—HO,+ AA,— BAgz 
My, = (H—pF)6+(B—pA)u—H®,+ FO,— BA,+AA,z 


(26) 


where 

A = [kg k(—y*)+k, Bk(py+yA—pA—)VA, 

B = [Kt kg kily—?) +k, I kgly-+-py—p—VA, 

C = [ky Q(.—p*) +1 ke 2yp—y*—1) V/A, 

D = [ky k§A—p?A) + kyl? p+p—2y) V/A, 

F = [kg k,(1—A*) +k, kg (2yA—y*— )/A, 

H = [k4ky(p—d%p) +k, BRA+yA—2y) V/A, 

A = k{(1—p?—A®+A%p*) + (1 —2y*7+-y4)+ 

+h ki kq(p— 2y?+ 4ry—2Ap+ 3py—2y*pA—2). 

4.6. The chords in this case will be subjected to different axial loads 

P, and P, which are assumed to vary sinusoidally. The differential 


equations (13) must in this case be written down for each chord, giving 
four simultaneous differential equations. The solution of these equations 
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are obtained in the same way as for the case treated above, and will be 
as follows: 


a,[} 2 [u*rnt-+3CL—2wP, B,]-+2wP, ¥ a,Bye+ 
+b,[4A L—2wP, YoBn]+ 2wP, Yo 2 b.Bns 
= {DL ya,+4BLyb, + J {tal +Ce)—0p(B-+De) + 


1 
+ k [A(uy+ywg)— B(wgtyu,)]+ 


+ 5 [Ola + Mop) —Diwoy +r} — 
—48Q[(A + B)/e+C+ D]/wn5 
a,[}4L—20P, YoBq]+20B, Yo S 4,Byat 
+b,[}ET tan! +4GKwnn? + $F L—2uP, r$B,]+ 


20 P16 d O.Bns 


1 BLya, + 4H Lyb, + = |0,(P —Ae)—0g(H + Be) + 
Ww 





1 
+ , [ F(uy+-ywp)—H (wg+yu,)|+ 


1 
+ E, [A(w,+Aw,) — Bw +Auy)) = 


—48Q[(F+-H)/e+A+-C]/wind 
pa,[ $2 Iu*nnt+4CL—2wPpBq]+2wPy pS a, Bny+ 


+ pb, [$4 L—2wP pg Yo By |+2wPx You p3 bs Bns 


- jDLa,+4BLb, +—|—04(B-+De) +094 +Ce)— 
Ww 


1 
— |, [Beer t+y0p)—A(wpt+yy)]— 


3 


1 
— Z [Dita +Awg)— Own +d) + 
1 





+48Q[(A + B)/e+C+D)/win 
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where p = k,/k, and k, = moment acting at A to cause unit rotation at 
A and k, = moment acting at A to cause unit rotation at B. 











4.5. The following assumptions are now made: 
Up= py = pu; dy = phy = pd. (25) 
The correctness of these assumptions will be discussed later. Using (23) 
and (25) the reactions due to the elastic medium can be expressed by 
means of the solution of (24), as: 


R, = —(A—pB)¢—(C—pD)u+A®,— BO,+CA,—DA, 
Ry = (B—pA)6+(D—pC)u— BO,4+-A®,—DA,+ CA, : 


(26) 
= —(F—pH)$—(A—pB)u+ F®,—H®,+AA,— BA, | 
My = (H—pF)$+(B—pA)u—HO,+ FO,—BA,+AAz 


where 

A= [ki ks ki 1—y*) +k, Kg ky(py+yA—pA—)/A, 

B= [kik ki(y—y) +h, Bkyly+pry—p—a)/A, 

C = [ky k§(1—p*) +47 k,(2yp—y*—1)/A, 

D = [ky K(A—p?A) + i kealy?p+-p—2y)/A, 

F = [kg k,(1—2*) + ky IB K2A—y2*— YA, 

H = [kg k,(p—A*p) +k, Wg A+ y°A—2y)]/A, 

A = Ki(1—p?—A?-+A%p*) + RG — 2? +4) + 

+k, kj k4(p— 2y?+ 4Ay— 2Ap+ Bpy—2y2pA— 2). 

4.6. The chords in this case will be subjected to different axial loa 

P, and P, which are assumed to vary sinusoidally. The differentia! 


equations (13) must in this case be written down for each chord, giving 
four simultaneous differential equations. The solution of these equations 
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are obtained in the same way as for the case treated above, and will be 
as follows: 


a,{} £Ia*n*+4C0L—2wP, B,]+2wP, ¥ a,Bqs+ 


+b,[44 L—2wP, YoBn|+ 2wP, Yo 2 b.Bns 


i ion mat ne 


9 
{DLpa, + 4BLpb, + =|04(A-+Ce)—0q(B+De)+ 


1 
+ 5 [Aust 109) — Bop ty) + 


Apres 3 tet hed ee 
nhac ae ee ee 


Soh 


+ ¢[O(m+dug)—Diwey+)]} — 


—48Q[(A+ B)/e+C+D)/w'n® 


a,[4AL—2w A YoBy|+ 20k, Yo > a,B,,+ 


Ney 
ae ney 


Sk it a ed ig ER. aa 
“ hee g 


+6,[4 ET wan +4GKwrn?+4F L—QwP, 2 B,.)+ 
+2wP, r > b.Bns 


$B Lpay-+ 4H Lp, +—|04(F—Ae)—Oy(H + Be) + 
] 
2 [ F(w, +ywp)—H(wpt+yuy,)|+ 
3 


+f [Ale + Arg) — Berg +duy)]} — 
1 


—480[(F +H)/e+A+C}/aens 
pat, | 4 E Iw*rn*+- $CL—2w BB n|+2wPy > @,B,,+ 
+pb, [4A L—2wP pg yo 8, |+2w BYoh > b.Bns 


— }DLa,+4BLb, +—~-|—6,(B+De)+6,(A-+Ce)— 
¥ wn B 


1 
re [ B(uy,+-ywg)—A(wp+yu,)|— 
3 


mm E [D(w,+Awg)— Own +d} + 


Medes ae tay eee tae ee SERRA Ra EDS gt rey - eS 
Ss Spal path ae deg se Sie act bat ENON 5 Se ee ™ 4 





+480[(A+ B)/e4+C+D)/win® 
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p4,[4AL—2wPz yo 8, |+-20Pz you 2 a, Bast 
+yb,[}ETu*ant+4GKwnn?+4F L—2wP5r28,]+ 
+-2wPp rip > b. Bus 
= }BLa,+4HLb,+—|—0,(H+ Be)+ (P+ Ae)— 


l 
— 5 Aut e)— Pept) 


] 
— 7 [Bley +p) —A wn} + 
1 


+-480[(F + H)/e+A+ B)/w'n’ 
(a == 1,3, 5,...; ¢ = 1,3, 5,...; » #8), 


4.7. For the case when the chords have negligible torsional and 
warping rigidities, the solution (27) simplifies to: 


a,| 4B Iw*an + tek, L—2wP, B,|+20P, ¥ O,Bns 





= fedk, Lua, 4-2ek, Ay ,/wn—48e'k, O/w'n> + 2w,/wn 
a,[ 4B Too*nn' + ek, L—2aPyBq]+ 20 Py S 4, Bug 


= }eAk, La, + 2ek, Ayp/wn+ 48e'k, Q/w'n> + 2wp/wn 
(n = 1,3,5,...; ¢ = 1,3,5,...; n 4 8) } 
where ¢ = 1/(1—A?); e’ = 1/(1—A) 
Aya = Sgr—AAzpi: Aye = Apr —AAgr- 





5. Convergency of the Fourier series 

It is not possible to give a general rule about the numbers of terms of 
the series which must be used in order to obtain results of sufficient 
accuracy for design purposes. Many factors influence the convergence of 
the series, and calculations will in each case clearly show the effects on 
the results of the number of terms in the series used. The following 
discussion is, therefore, intended to illustrate how much work is involved 
in using the formulae of this paper in actual design. 


5.1. One of the most important factors to influence the convergence 
of the series is the value of the axial load in a chord compared with the 
corresponding buckling load for the straight chord. Table 1 shows clearly 
how the accuracy of the results is affected by the magnitude of the axial 
load. For small axial loads, the maximum deflexions are overestimated 
by using an odd number of terms in the series. When, however, the axial 
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TABLE 1 





Percentage error in 


Axial load as maximum deflexion 


percentage of 
buckling load One term Two terms 








° +744 —o'79 

9°5 +712 —o"6g 
19 +661 — 0-60 
38 +473 — O44 
47 +313 — 0°39 
57 +0°83 —0o"40 
76 —7°53 — 1-00 
85 —15°59 — 2°33 











loads are very large, the even terms in the series will change sign, and 
the deflexions will be underestimated when an odd number of terms are 
used. For the particular case shown in Table 1, this takes place when 
the axial load is about 60 per cent of the corresponding buckling load 
for the straight chord. 


5.2. Another important factor which influences the convergence of the 
series is the magnitude of the flexural rigidities of the chords compared 
with the stiffness of the elastic medium. It can be seen from equation (17) 
that when £1 is large compared with k, the convergence of the coefficients 
will be rapid. Table 2 shows how the convergence of the series for u 


ay 


TABLE 2 





Percentage errors in the maximum deflexions and rotations 





a GK large, EY small GK small, ET large GK and ET small 


terms u“ ry “ od u ¢ 


+ 10°43 —o18 +445 —O°35 + 4°02 — 2°02 
I'l4 +0O°72 —O'54 +o°387 o’31 +1°74 

+ 0°23 -0°36 +O'12 —0O'22 +-0°06 — 318 
0°05 +o'10 — 0°03 +O°05 —O'or +0°54 

+ 0°03 —O13 + 0°02 — 0°03 +0°Or — 2°38 





























and ¢ is affected by the magnitudes of the torsional and warping 
rigidities of the chords compared with the stiffness of the elastic medium, 
The axial load in this case is about 30 per cent of the corresponding 
buckling load for the straight chord, and the flexural rigidity the same in 
all cases. Results obtained by using six terms are assumed to be exact. 


5.3. For design purposes, the moments and torques must be determined 
from the deflexion and the rotation curves. This can be done either by 
the differentiation of the series for these curves, or by considering the 


5092 .53 H 
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equilibrium of the chords, i.e. by the integration of the series. The latter 
method will, generally, give the more accurate results. Table 3 indicates 
the zceuracy of the two methods. The bending moments at midspan are 
here calculated for a chord when the effect of torsion is negligible, and 
the axial load is about 40 per cent of the corresponding buckling load for 
the straight chord. 


TABLE 3 





Percentage error in moment at midspan 





By differentiation By tntegration 





+45°43 +8-26 
— 9°86 — 2°32 
+3°81 +0°49 
— 1°03 — 060 
+0°96 +0°07 
—0O°44 ° 











5.4. In practice, the maximum stresses in the chords will usually occur 
when the cross-girders are symmetrically loaded. When torsion and warp- 
ing are negligible, the bending moments can usually be calculated with 
sufficient accuracy by using three terms in the series. This means that 
the evaluation of the results will only consist in the solution of three linear 
simultaneous equations. For the cases when torsion and warping are 
taken into account, the convergence may sometimes be rather slow, 
especially if an accurate estimation of the torque is required. For cases 
where the torque plays an important part, however, the values of the 
torsional and/or the warping rigidities will be large, and the convergence 
consequently rapid. Usually three terms in the series will be sufficient, 
and the calculation of the results will consist in the solution of six linear 
simultaneous equations. This can be done on an ordinary desk calculating 
machine. 


6. Application of the theoretical analysis to actual design 

It is concluded from the above discussion that the theoretical results 
will, for most practical cases, provide the necessary data for the design 
of compression chords for through-bridges without too many computa- 
tion being involved. Before applying the theoretical results to the design 
of actual bridges it will, however, be of importance to consider to what 
extent the assumptions on which the analysis is based are fulfilled in 
actual bridge-structures. 


6.1. The theoretical analysis assumes linear elastic behaviour of the 
structures. In welded structures this may not be far from being true, but 
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with bolted or riveted connexions this assumption will not be completely 
fulfilled. It has been shown, however, by experiments with a riveted 
aluminium model (3), that good agreement was obtained between the 
experimental and the theoretical results. This indicates that, even for 
bolted or riveted connexions, the theory will predict the behaviour of 
the chords accurately. It is, however, of great importance to base the 
theoretical calculations on reliable data, especially with regard to the 
behaviour of the U-frames. The assumption that only the elastic behaviour 
of the chords and U-frames is considered will probably cover most practical 
cases since it is generally the aim of designers of bridges to keep the 
maximum stresses well below the yield point in order to prevent any 
danger of fatigue. 


6.2. The compression chords are assumed to Le perfectly straight before 
any loading is applied to the bridge. This will generally not be the 
case in practice. In studying the buckling behaviour of the straight chords 
any initial departures from straightness will lead to serious reductions of 
the buckling loads. When deflexions and rotations are considered for a 
loading applied to the cross-girders, the initial imperfections will, however, 
usually be very small compared with the deflexions and rotations produced 
in the chords by the deflexions of the cross-girders. In addition, the axial 


loads in the chords will be small compared with the corresponding 
buckling loads for the straight chords. The neglect of initial departure 
from the straightness of the chords, therefore, in most practical cases will 
not affect the accuracy of the results appreciably. 


6.3. The whole of the theoretical study is based on the assumption that 
the lateral stiffness of the U-frames can be replaced by an equivalent 
elastic medium. For this to be true in practice, the U-frames must be of 
equal stiffness and be evenly spaced throughout the span, and the number 
of U-frames must not be too small. Calculations show, however, that 
with only four intermediate U-frames the results calculated from the 
assumption of a continuous elastic medium are very nearly equal to those 
calculated from the assumption of evenly spaced discrete elastic supports. 
The number of U-frames will usually in practice be quite large, and the 
assumption of a continuous elastic medium can safely be made. 


6.4. For the cases of unsymmetrical loading on the cross-girders it was 
assumed in section 4.5 that up, = wu, and ¢, = pd,. In very extreme 
cases, e.g. for bridges having very flexible main girders subjected to high 
wind loadings and with cross-girders subjected to very high unsym- 
metrical loadings, these assumptions will not be correct. Experiments 
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with model bridges (3) have shown, however, that even for loadings which 
go far beyond the conditions which are likely to be encountered in practice, 
the assumptions lead to results which in most cases should be accurate 
enough for design purposes. 


6.5. The behaviour of the chords has been studied both for the cases 
when the effect of torsion is neglected and when it is taken into account. 
It is difficult to give a general rule as to when the effect of torsion 
becomes so important that it should be taken into account, but the effect 
of neglecting it will be as follows. In the case when the chords are of 
closed box section, the high torsional rigidity will reduce the deflexions 
and rotations of the chords as calculated on the assumption that torsion 
can be neglected; and will consequently lead to too safe designs. In the 
case when the chords are of open section, the high warping rigidity will 
increase the deflexions and rotations of the chords calculated on the assump- 
tion that torsion can be neglected, and will consequently lead to unsafe 


designs. 


6.6. It can be concluded from the above discussion that the analysis 
presented in this paper should, within the limits set by the assumptions 
which have been made, be applicable to the design of most of the usual 
types of through-bridges. 
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SUMMARY 
The wave functions in some axisymmetric scalar diffraction problems involving 
one or more circular disks are obtained as contour integrals to show that the 


solutions of these problems depend on the solutions of Fredholm integral equations 
of the second kind in one independent variable. 


1. Introduction 


In 1956 D. 8. Jones (1) developed a method of reducing to the solution 


of a one-dimensional Fredholm integral equation of the second kind the 
problem of the diffraction of harmonic sound waves by a perfectly rigid 
circular disk, the incident waves being symmetric about the axis of the 
disk. In this method the perturbation wave function, that is, the difference 
between the total and incident wave functions, is regarded as being due 
to a fictitious distribution of wave doublets over the disk, the boundary 
condition on the disk leading to the Fredholm equation governing the 
problem. Work by Green and Zerna (2) and the author (3, 4, 5, 6) on the 
solution of axisymmetric potential problems for circular disks and spherical 
caps (using contour integral representations of potential functions to 
reduce such problems to the solution of integral equations) suggested that 
this approach might be applied to the diffraction problem for the perfectly 
rigid disk to give a simpler derivation of Jones’s integral equation; and 
on investigation this proved to be so. Further, it was found possible to 
apply the contour integral method to the diffraction problem for a per- 
fectly soft circular disk and to diffraction problems for two or more disks. 

It was, however, pointed out by a referee of the original version of this 
paper that the basic technique and the results for diffraction problems 
for a single disk by the contour integral method had already been antici- 
pated by Bazer and Brown (7) in a recent New York University report, 

(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961) 
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in which they consider axisymmetric diffraction problems for a screen 
containing a circular aperture. By Babinet’s principle the results for these 
problems can be applied to diffraction problems for a single disk. 

Bazer and Brown, however, take as their starting-point contour integral 
representations of the perturbation wave functions without indicating 
how these representations can be derived from more elementary considera- 
tions. This defect in their approach can be overcome and in sections 2 
and 3 of this paper it is shown that, by regarding the perturbation wave 
functions as being due to fictitious distributions of wave sources and 
doublets over the disk, we can derive the required integral representations 
by means of a solution of the wave equation given by Bateman (8). 

The application of the contour integral method in diffraction theory is, 
however, not limited to problems for a single circular disk; it is possible 
to apply the method to problems for two or more parallel disks, or screens 
containing circular apertures, provided that the incident field is axisym- 
metric. In section 4 the diffraction of harmonic waves by two perfectly 
rigid parallel circular disks is considered and the governing Fredholm 
equations of the problem found. For incident plane waves whose wave 
number is small an approximation to the scattering cross-section of the 
disks is obtained. Results for the diffraction of harmonic waves by two 
perfectly soft parallel circular disks are stated in section 5. In a subse- 


quent paper the diffraction of harmonic waves by two perfectly rigid 
screens, each containing a circular aperture, will be considered. 


2. Representations of the Helmholtz functions for rings of sources 
and doublets 

In the problem of the diffraction of harmonic waves by one or more 
perfectly rigid circular disks we regard the perturbation wave function, 
that is, the difference between the total and incident wave functions, as 
being due to fictitious distributions of wave doublets over the disks. If the 
disks have a common axis about which the incident field is symmetric, 
the densities of the doublet distributions are symmetric about this axis, 
the axes of the doublets being parallel to it. Similarly, in the problem 
of the diffraction of harmonic waves by a system of perfectly soft circular 
disks having a common axis about which the incident field is symmetric, 
the perturbation wave function is regarded as being due to fictitious axi- 
symmetric distributions of wave sources over the disks. 

We therefore begin by finding an integral representation for the wave 
function due to a ring of sources. We make use of an expression given by 
Bateman (8) for a wave function which is symmetric about an axis and 
finite on the axis. Using cylindrical polar coordinates (a, 4,2) referred to 
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some origin O and the axis of symmetry as z-axis, the wave function 
O(a, z,t), which is symmetric about the axis w = 0, and has the values 
P(z,t) at points on it at time t, is given by 


27 
O(a, z,t) = x ( P(:-+ie cosy, t—“ sing) dy, (2.1) 
‘ 


where c is the wave velocity. For harmonic waves we write 
O(a, z,t) = V(w, ze, 


where w is the frequency of the waves, so that V(w,z) is a solution of 
Helmholtz’s equation 


a2 2 

ee P adtad ve = +i8)7 (w,2) = 
022 Oo”  wlo 

k = w/e being the wave number. Writing P(z,t) = p(z)e, we have from 


(2.1) the Helmholtz function V(w,z), which has the values p(z) on w = 0, 
and is given by 
2r 


a | p(e-+ier 008 jexp(—iker sin) dip 
7 


0 


1 
7 


[ me +i cos s)cos(kar sin ys) dis 
) 


. 
( 


+ pers ro eae p(z+it) dt, (2.3) 
w 

where the final integral is obtained by use of the substitution t = wcosy. 

For a ring of wave sources of radius pin the plane z = ¢, the Helmholtz 


function v,(@, z; p, ¢) at the point (a, ¢d, z) is 


d¢’, 
(2.4) 
the total source strength of the ring being unity. On the axis w = 0 
exp[ —ik{(p?+(z—{)"}'], 
ip (2—o ; 





- = {ex xp Gn... 24 52 2earp cos(d—d’) + (z—C)*}4] 
v(m, 2; p,f) = + p?— esp cos(d—¢') +(z —f)t 


2a 





v,(0,2; p,f) = plz) = 


hence from (2.3) we obtain 





f cos k(w*—t)texp[ —ik{p 2+ (z—{+it)*}] dt. 


1 
vile, 25 pt) = = (o?—#)Hp eer 


. 
—-w 
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From (2.4) we have 


vu, (a, z; p> t) — 0, (p, <e w, 2), 
and hence obtain as a second integral representation for v,(w,z; p, ¢) the 
expression 
cos k(p?—t*)* exp[ —ik{a?+ (z— [+ it)?}4] dt. 
(p?—t?)? + (z—C+-2t)?}# 
—P 


We next consider a ring of wave sources of total strength —q in the 
plane z = { and a ring of sources of total strength q in the plane z = (+-8{, 
both being of equal radius p. The Helmholtz functions for these rings are 
then obtained from (2.5). In the limiting case when g > o and df > 0 in 
such a way that ¢3{ — 1, the Helmholtz function v,(m, z; p, ¢) at the point 
(a, }, 2) due to a ring of doublets of radius p in the plane z = ¢ is given by 


p 
i __ & f cosk(p?—t*)! é fexp[ —ik{w?+ (z—-+ it)?}*] s 
neti = 2 | car al erect |e Oo 





1 
v,(a,2z; p,f) = oe 


(2.5) 





om 
the total strength of the ring being unity and the axes of the doublets being 
parallel to the z-axis. 

Equations (2.5) and (2.6) are the bricks out of which the Helmholtz 
functions in diffraction problems for circular disks are to be built. 


3. The diffraction of harmonic waves by a single disk 

We now consider the diffraction of harmonic waves by a circular disk 
of radius a and take the centre O of the disk to be the origin of cylindrical 
polar coordinates (a, ¢, z), the axis of the disk being taken as z-axis. The 
disk is thus given by z = 0 (0 < w <a). The Helmholtz function V,(a, z) 
for the incident field is assumed independent of ¢; hence the total Helm- 
holtz function V,(a@,z) is also independent of ¢ and is given by 

Vi(w, z) = Vi(w, z)+V (oe, z), (3.1) 

where V(a,z) is the perturbation Helmholtz function. This latter func- 
tion must be continuously differentiable at all points except possibly 
those on the rim of the disk, where it is finite. At large distances from the 
disk V(w,z) satisfies the Sommerfeld radiation condition (9), that is, it 
represents diverging waves. If the disk is perfectly soft, the total Helm- 
holtz function vanishes on the disk; hence 


V(w,z) = —Vi(w,z) = e(w), say,onz=0(0<am <a). (3.2) 

If the disk is perfectly rigid, the normal derivative of the total wave 
function vanishes on the disk; hence 

eV(m,z) __ Mla, 2) a 


‘ R, =0(0<wca). (3.3 
= —— =flw), say,onz=0(0<w <a). (3.3) 
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The functions e(w) and f(a) are supposed continuously differentiable even 
functions of w. These conditions are sufficient to characterize the func- 
tion V(w,z) uniquely (10). 

For a perfectly soft disk the perturbation function V (a, z) is regarded as 
being due to a fictitious distribution of sources on the disk, the distribution 
being symmetric about the centre of the disk, so that its density can be 
taken as o(p) per unit area at the point (p, ¢’, 0) of the disk. We thus have 


V(w,z) = { 2npa(p)v, (a, z; p, 0) dp, 
0 
which, by using (2.5) and interchanging the order of integration, gives us 
iw) f (3.4) 
hey age R . 


where R = {w?+-(z+1t)*}! 


and g(t) is a complex continuous even function given by 


a 


gt) = 4 | 


po(p)cos k(p? —t?)* 


(pt?) dp (t > 0), 





=g(-—t) (t< 0). 


Similarly for a perfectly rigid disk the perturbation function V(m, z) is 
regarded as being due to a distribution of doublets over the disk, the axes 
of the doublets being parallel to the z-axis. Since the distribution is sym- 
metric about the centre of the disk, its density can be taken as 7(p) per 
unit area at the point (p, 4’, 0) of the disk. We thus have 

. 
V(w,z) = | 2npt(p)v,(or, z; p, 0) dp. 
0 
Using (2.6), interchanging the order of integration and carrying out an 
integration by parts, we find that 


a 
1 
2% 
—6 


Vi(w,z) = (3.5) 


where j(t) is a complex continuous odd function given by 


: d  pr(p)cos k(p?—t2)! 
t) == 4— 
jt) = 4 (eh dp 
t 


(t > 0), 





—j(—t) (t< 9). 
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From (2.4) we have 
v,(w, a p> f) ae v,(p, e w, z), 
and hence obtain as a second integral representation for v,(w,z; p,¢) the 


expression 
Se ate 2 ve s#)2 
v,(w, 2; p,f) = 4 ea etre i+) dt. (2.5) 
—p 





(p?— 0?) ow? + (z—£+ it) 


We next consider a ring of wave sources of total strength —gq in the 
plane z = { and a ring of sources of total strength qg in the plane z = (+-8£, 
both being of equal radius p. The Helmholtz functions for these rings are 
then obtained from (2.5). In the limiting case when g > o and 8f > 0 in 
such a way that ¢8{ > 1, the Helmholtz function v,(m, z; p, f) at the point 
(ar, }, z) due to a ring of doublets of radius p in the plane z = { is given by 


p 
; ‘ __ & f cos k(p?—#*)! @ [exp[—ik{w*?+ (z—+-it)?}4] 
v,(@, 2; p,¢) = — | (p2—t?)! | {a?+-(z—{-+-it)?}! 


7 
the total strength of the ring being unity and the axes of the doublets being 
parallel to the z-axis. 

Equations (2.5) and (2.6) are the bricks out of which the Helmholtz 
functions in diffraction problems for circular disks are to be built. 





dt, (2.6) 


3. The diffraction of harmonic waves by a single disk 

We now consider the diffraction of harmonic waves by a circular disk 
of radius a and take the centre O of the disk to be the origin of cylindrical 
polar coordinates (a, ¢,z), the axis of the disk being taken as z-axis. The 
disk is thus given by z = 0(0 < w <a). The Helmholtz function V,(z, z) 
for the incident field is assumed independent of ¢; hence the total Helm- 
holtz function V,(@, z) is also independent of ¢ and is given by 

Vi(@w, z) = Vi(w,z)+V (a, 2), (3.1) 

where V(m,z) is the perturbation Helmholtz function. This latter func- 
tion must be continuously differentiable at all points except possibly 
those on the rim of the disk, where it is finite. At large distances from the 
disk V(mw,z) satisfies the Sommerfeld radiation condition (9), that is, it 
represents diverging waves. If the disk is perfectly soft, the total Helm- 
holtz function vanishes on the disk; hence 


V(w,z) = —Vi(aw,z) = e(w), say,onz=0(0<aw<a). (3.2) 

If the disk is perfectly rigid, the normal derivative of the total wave 
function vanishes on the disk; hence 
eV (w, z) = _ Vole, 2) = f(a), 


say,onz=0(0<aw<a). (3.3 
a Op y (0<@ <a). (3.3) 
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The functions e(w) and f(a) are supposed continuously differentiable even 
functions of w. These conditions are sufficient to characterize the func- 
tion V(w,z) uniquely (10). 

For a perfectly soft disk the perturbation function V (a, z) is regarded as 
being due to a fictitious distribution of sources on the disk, the distribution 
being symmetric about the centre of the disk, so that its density can be 
taken as o(p) per unit area at the point (p, ¢’, 0) of the disk. We thus have 


a 


V(w,z) = | 2npo(p)v,(ar, 2; p, 0) dp, 
0 


which, by using (2.5) and interchanging the order of integration, gives us 
a 


V(w,2z) = : | = dt, (3.4) 


— 


—a 
where R = {w?+(z+1t)*}! 
and g(t) is a complex continuous even function given by 


g(t) = 4 ( teen OY ap 


(ee sat 
t 





=g(—t) (t< 0). 


Similarly for a perfectly rigid disk the perturbation function V(a, z) is 
regarded as being due to a distribution of doublets over the disk, the axes 
of the doublets being parallel to the z-axis. Since the distribution is sym- 
metric ahout the centre of the disk, its density can be taken as r(p) per 
unit area at the point (p,¢’,0) of the disk. We thus have 


V(w,z) = | 2mpr(p)v,(a, z; p, 0) dp. 
0 


Using (2.6), interchanging the order of integration and carrying out an 
integration by parts, we find that 


a 


V(w,z) = o | = dt, (3.5) 


—a 


where j(t) is a complex continuous odd function given by 


t 





(p?—t*)! 


—j(—t) (t< 0). 
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The square roots in (3.4) and (3.5) are interpreted as follows. For z > 0 
we have R= R,e¥, (3.7) 


and for z < 0 = f (3.8) 
where R, : 0<t<a, 
< 


< 
-6<0 for—a<t<0. 
When z = 0, we have 
R= (w*—?)' (w > |t)), 

while 

R= +i(?@—am*)! (wo < |t|) ifz+0 through positive values, <o 

= Fi(®—m’)! (w < |t|) ifz+0 through negative values, — 

the upper signs holding for 0 <¢ <a and the lower for —a <t < 0. 
There is thus a discontinuity in R at z = 0 (w < |t|). The integral in 
(3.4) is to be interpreted as a Cauchy integral at the point (z = 0, a = 0). 

Equations (3.4) and (3.5) give the integral representations used by 
Bazer and Brown as their starting-point. In either case V(w,z) satisfies 
the Sommerfeld radiation condition and, provided g(t) or j(t) is a con- 
tinuous function, can be shown to be a continuously differentiable func- 
tion satisfying Helmholtz’s equation (2.2) at all points except those on the 
rim of the disk where it is finite. The condition (3.2) or (3.3) satisfied by 
V(a@,z) on the disk leads to a Fredholm integral equation of the second 
kind satisfied by g(t) or j(t), the derivation of this being given by Bazer 
and Brown (7). Since the method of deriving these integral equations is 
to be extended to problems for two disks, an outline of the derivation of 
the equation satisfied by j(t) is now given. 

If the disk is perfectly rigid, we require V(a, z) to satisfy the condition 
(3.3) on the disk. Since 


we have 
éfexp(—ikR)\ _ ae eS exp(—ikR))\, 
oe) ~ et 


so at any point not on the disk we have 


iw 





eV (a, 2) * it)(2 +-it)exp(—tk R) di 
2 tél R ‘ 


whence 





2i [ pTUe=) dp — fi (n| SERRE) expt Fik(e +i} dt, (3.10) 


R 


C2 
0 “a 
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the upper signs in the integral on the right-hand side holding for z > 0 
and the lower for z < 0. It can be shown that the limit of this integral 
as (w,z) approaches a point on the disk is equal to the value of the 
integral at that point. Hence using (3.3) and (3.9) we find that on the 
disk z= 0(0 < w <a) 


h(a) = { ef (p) dp 
0 


dt + 








ah f eos bw — a; { Gear ='y 


(or? —t?)* (t?—m?)! 


« 


0 
a 


+i { j(t)sinh kt dt. (3.11) 
0 
We now make use of the following result due to Jones (1). The solution 
of the integral equation 


w 


(ae 2__ ¢2)t dt=diw) (0<w<a), (3.12) 


(wi 





a d(ar)cosh k(t? —aw*)* 


(t?—?)t da (0 < t< a). (3.13) 





If we therefore regard (3.11) as a Volterra equation for j(t) in which the 
first term on the right-hand side is equal to a known function, we find 
from (3.13), after some manipulation, that its solution is 


j(t) 





n(t) (—a <t <a), (3.14) 


Bt i [ j(8)sinh k(t—s) ae 
= cnn se = 


T 
. 
—a 


where n(¢) is a continuous odd function given by 


t 
2d a 


aa d 0<t<a 
~ it di 2—w)t es 





= —n(—t) (—a <t< 0). 


Equation (3.14) is a Fredholm integral equation of the second kind for 
j(t), n(t) being a known function. From the theory of Fredholm equations 
(11) the eigenvalues of (3.14) are real and hence the equation possesses 
a unique solution. We may note that the eigenfunctions of the kernel of 
(3.14) are the angular spheroidal wave functions of order zero (10). The 
equation for j(t) given by Bazer and Brown is rather more complicated 
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than (3.14), but reduces to it after some manipulation. Equation (3.14), 
however, has the same kernel as the equation obtained by Jones (1), 
although the unknown function j(¢) in (3.14) is not to be identified with 
the unknown function in Jones’s equation. 

Various methods are available for solving (3.14). When k is small, the 
most convenient method appears to be that of iteration, and in this way 
Bazer and Brown obtain an approximate solution when the incident 
Helmholtz function represents a normally incident plane wave. In the 
iterative solution of (3.14) the first approximation is obtained by putting 
j(t) = n(t), that is, the integral on the left-hand side of the equation is 
neglected. This approximation can be shown to have a direct physical 
significance by considering the non-radiating problem for the disk. Here 
we assume there is no radiation from the disk and thus replace the Sommer- 
feld condition on V by the condition that V represents standing waves at 
infinity, the remaining conditions being as before. The form for V is then 


a 


2iV(w,2z) = [ 


—a 


jo(thcos kR 


R ™ 


where j,(t) is a continuous odd function, and, proceeding as before, we 
have on the disk 


w 7 
” 


| f(p) dp = hw) (0< aw <a). 





tj,(t)cos k(a? —t?)* Page 
(ar? — t?)# = 
. 0 


On solving this equation we find that j,(¢) = n(t), so that j,(t) is identical 
with the first approximation to j(f) obtained from the integral equation 
(3.14). 

If the disk is perfectly soft, V(w,z) satisfies the condition (3.2) on the 
disk. By a method similar to that used for a perfectly rigid disk the 
integral equation satisfied by g(t) is found to be 


a 


ait)—* [ g(s)sinh k(t—s) ds = m(t) (—a <t <a), 
7 t—s 


when m/(t) is a continuous even function given by 


t 
_ 2d f we(w)cosh k(t?—w)! 
oa (?@—o)t 
0 





(0 <t <a) 


=m(—t) (—a <t < 0). 
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It follows from the theory of Fredholm equations that equation (3.15) 
uniquely determines a continuous even function g(t). 


4. The diffraction of harmonic waves by two perfectly rigid parallel 
circular disks 


The method of the previous sections can also be applied to the diffrac- 
tion ef harmonic waves by two perfectly rigid parallel circular disks, a 
problem which does not so far seem to have been considered. We suppose 
that two disks &, and &,, of radii a and b, are parallel and at a distance 2f 
apart, the line joining their centres being perpendicular to the planes of 
the disks and therefore the axis of symmetry for the configuration. 
We take this line as z-axis for cylindrical polar coordinates (w,¢,z), the 
origin O being at the midpoint of the line. In addition we use cylindrical 
polar coordinates (w,¢,z,) and (w,¢,2,) with the centres O, and O, of 
the disks 2, and &, respectively as origins, so that z, = z+-f and z, = z—f. 
The equations of the disks £, and &, are thus given b, . = —/f, z, = 0 
(0 <w <a) andz=f,z, = 0(0 <m < dD) respectively. 

We suppose that the incident Helmholtz function V,(a, z) is independent 
of ¢, so that the total Helmholtz function V,(w, z) given by 


V;(ar, 2) = Voor, z)+V (er, z) (4.1) 


is also independent of ¢. The perturbation Helmholtz function V(m, z) is 
continuously differentiable at all points except those on the rims z = —f, 
@o—a,andz=f, w = b, where it is finite; it satisfies the Sommerfeld 
condition at large distances from the disks. Since the disks are perfectly 
rigid, the normal derivative of the total Helmholtz function vanishes on 
either disk. We thus obtain 


(0 <a <a), 


= 0 (0<mw<b). (4.2) 


The functions f,(@) and f,(a) are taken to be continuously differentiable 
even functions of w. 

The perturbation function V(w, z) is regarded as being due to fictitious 
axisymmetric distributions of doublets on the disks, the axes of the 
doublets being parallel to the z-axis. To obtain an integral representation 
of V(w,z) we therefore add together the Helmholtz functions for the two 
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distributions and by a method similar to that by which (3.5) is derived 
we find that 
a b 
2iV (w,z) = sn \ ees ae dt, (4.3) 
R, R, 
=b 








2 
—a 


where R, = (w*+(z,+it)}', By = (ew? + (z,+it)*}, 


and j,(t) and j,(t) are complex continuous odd functions of ¢ related to the 
densities 7,(a7) and 7,(a) of the doublet distributions on L, and L, by 
expressions similar to (3.6). The square roots R, and R, are defined by 
(3.7), (3.8), and (3.9) with z replaced by z, and z, respectively and a 
replacing 6 in R,. 

The proof that V(a,z) satisfies the required conditions other than (4.2) 
is similar to that for the corresponding problem of a single disk. From 
the conditions (4.2) which are satisfied by V(w,z) on the disks we now 
derive two Fredholm integral equations of the second kind for j,(¢) and 
),(t). Since ‘ 
ay 1 it) Bn (n = 1,2), 

Ow 


it follows that at any point (w,z) not on either disk 


2ia —. 
Oz 
a b 


ss | f Salles Spee) at f 2afeart ert a, 








ow 1 J 2 
—a —b 


so that 


ai [ pea = | jyn| Se erp NY Fexpt Fike, +i0} dt + 
ny 1 


Cz 





. 
a 


b 


+ | Jol | Seep Ane) exp Fik(=, +i} dt. (4.4) 
2 


. 


The upper signs in the first integral on the right-hand side of this equation 
hold for z, > 0 and the lower for z, < 0, while in the second integral the 
upper signs hold for z, > 0 and the lower for z, < 0. 

We first consider the condition on the disk 2, given by z = —f, z, = 0 
(0 <@ <a). As in the problem for the single disk the limit of the first 
integral on the right-hand side of (4.4) as (w,z) approaches a point on 
the disk =, is equal to the value of the integral at that point. We thus 
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find that for z = —f (0 < w <a) 
h,(w) = | pf (p) dp 
0 
s ’ tj, (t)cos k(w? —t?)! ditt tj,(t)sinh k(t?—aw*)! dt +. 
(w?—t?)* (t?—*)! 








+i | i(@sinh kt dt + 


0 





+ 
21 


b 
. {(—2f+it)exp[ —ikfa?+ (—2f+-it)?}*] 
J sx| fort (— prin J 
—b 
+exp{ik(—2f+ it)}| dt. (4.5) 


We note that, since z, = —2f on &,, the lower signs are to be taken in the 
last integral on the right-hand side of (4.4). 

We reduce the integral equation (4.5) to a Fredholm equation by regard- 
ing it as a Volterra equation for j,(t) in which the first term on the right- 
hand side is equal to a known function. After some manipulation we 
obtain its solution by means of (3.12) and (3.13) in the form 


if j,(8)sinh k(t—s) Pia 
t—s 





ji(t)— 


-_ 
a 
* 


—@ b 
he ase * jo(8)[ (t—s)sinh k(t —s) + 2if cosh k(t—s)] 
~ exp( 2ikf ) | Gata? ds 
—> 


=n,(t) (—-a<t<a), (4.6) 





where n,(t) is a known continuous odd function given by 





t 
2d E: k(t?—a?)! dw (0 <t <a), 


me) a 

0 

= —mn,(—t) (—a <t< 0). 

The condition on the disk £, given by z = f, z, = 0 (0 < w < b) gives 

a second Fredholm equation for j,(t) and j,(¢). This is found in the same 


way as (4.6) to be 


(@—w)} 


(4.7) 


b 
ja(t)—~ ( #e k(t—s) 4, 
Wy t—s 
—b 


_ tear) | Licata at eee ds 
7 








(t—s)*+-4f? 
=n,(t) (—b<t<b), (4.8) 
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where n,(t) is a known continuous odd function given by 


dw (0 <t<b), 





n,(t) = — — (P—w)t 


t 
2d [ wh,(w)cosh k(?—w"*)! 
nt dt 
0 


=-—n,(—t) (—b <t< 90), 


hy(w) | pf2(p) dp. 
0 


Equations (4.6) and (4.8) are Fredholm equations of the second kind 
for j,(t) and j,(t). From these equations j,(¢) and j,(t) are seen to be odd 
functions, whilst by a method similar to that used in the potential problem 
for two caps (5) it can be shown that j,(¢) and j,(¢) are continuous and 
uniquely determined by (4.6) and (4.8). Various methods can be used to 
obtain approximate solutions to the equations. For instance, when the 
disks are a large distance apart and the wave number & is small, the 
equations can be solved by iteration, the first approximations to j,(?) and 
jo(t) being n,(t) and n,(t) respectively. 

When the disks are of equal radius a, equations (4.6) and (4.8) are in 
the same basic interval (—a,a) and can be replaced by two equations, 
each determining a single function, in the following way. We write 


JL = HAMLIAO} 
with corresponding definitions for n,(t). On addition and subtraction of 
(4.6) and (4.8) we obtain 
a 
j.0-* | i.| 4 


t—s 
—a 


exp(— 2ikf ){(t—s)sinh k(t—s)+ 2if cosh k(t—s)} ds 
- (sap? | 


=n,(t) (-a<t<a). (4.10) 





Scattering of a normally incident plane wave by two equal disks 
We suppose the disks to be of equal radius a. The Helmholtz function 
for a normally incident plane wave travelling in the positive direction of 
2 is V(a@,z) = exp(—tkz), 
so that 
fi(w) = tkexp(tkf), f.(w) = ikexp(—ikf) (0 < mw <a). 


From (4.7) and (4.9) we find that 


9; SS 
n,(t) = * exp( ik )sinh kt, n,(t) = —exp(—ikf )sinh kt (—a <t <a). 
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We obtain the integral equations for the problem by putting these values 
of n,(t) and n,(t) into (4.6) and (4.8). It is, however, preferable to use 
equations (4.10) and for later convenience we give these in non-dimen- 
sional form. If we substitute 


t{=az, s=ay, «a=ak, f=a), 


and define new functions J,(x) by 
2; ¢ 
ji.) = — cos aAJ,(t/a), j()= Sale aA J_(t/a), (4.11) 
7 


we obtain the equations satisfied by J, (x) as 


a sinh a(x—y) 
“4(*) £ fu irae ay 





Fi Sen. y)sinh a(a—y)-+ 2iA cosh a(x— Y}) 
(yan | 


=sinhar (—l<2x<1). (4.12) 
The scattering cross-section of the disks is given by (12) as 
—(4n/k)im A, 


where A is the amplitude of the scattered wave in the direction of propa- 
gation of the incident wave. From (4.3) we find that at a large distance r 
from the origin O along the positive z-axis, V(w,z) has the form 


V(w,z) = = { [exp(—ikf)j,(t)+-exp(ikf )j2(t) sinh kt dt+- 
+O(r-%); 


hence the scattering cross-section S of the disks is given by 
a 


S = TP | [exp(—ikf)j,(t)+exp(ikf )j2(t) |sinh kt dt 
0 
16a? 


1 
me a wwe SR [ [cos*aA/, (2) +-sin*aAJ_(2x) |sinh ax dx. (4.13) 


a 


« 


0 
When the distance 2f between the disks is large compared with their 
radius, we obtain the solutions to (4.12) for incident waves of small wave 
number by iteration and, provided a and 1/A are sufficiently small, the 
power series so obtained are convergent. The first few terms in these 


5002.53 I 





1l4 

series are 

J,(2) = la) —* E(x) * exp(—2iad) (2) F exp(—2iad)Z,(z)— 
Se ar 


$exp(—6iad)L(2)+ ofa’ x 
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The scattering cross-section S of the disks is then found from (4.13) as 


S afl oh spl & 1 2 68 mt 
Sra® 219A? «SAA 27‘ 15\2 8lm\® 1575A4 


+f 324 44 OE J+ 
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5. The diffraction of harmonic waves by two perfectly soft parallel 
circular disks 

The diffraction of harmonic waves by two perfectly soft parallel circular 
disks can be treated in a similar way to the corresponding problem for 
two perfectly rigid disks. Only the results are stated here. 

We use the same coordinate system as in section 4 and take ),(a, z) 
and V(w,z) to be the incident and perturbation Helmholtz functions 
respectively, so that the total Helmholtz function V,(w, z) is given by 

Vi(@, z) == (ew, z)+V (a, z). 
We thus require a Helmholtz function V(a, z), which satisfies the Sommer- 
feld radiation condition at large distances from the disks and which is con- 
tinuously differentiable at all points except those on the rims of the disks 
where it is finite. Since the disks ©, and , are perfectly soft, we obtain 
V(w,z) = —Vi(w,z) on &, and X,, 

=¢(w~) onz=-—f,z,=0 (0<w <a), 

=e(w) onz=f,z.=0 (0<wmw < 5d), 
the functions e,(w) and e,(a) being continuously differentiable even func- 
tions of w. 

The function V(w, z) is obtained by adding two integral representations 
of the form (3.4) for distributions of sources on the disks, giving 


a 


b 
2V(w,2) = jee we (eee ate dt, 
1 
—b 








where the first and second integrals are to be interpreted as Cauchy 
integrals at the points (zc = —f, w = 0) and (z = f, w = 0) respectively, 
and #, and R, are as defined in section 4. Here g,(t) and g,(t) are complex 
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continuous even functions, uniquely determined by the Fredholm equa- 
tions 


,()—! jee aa 





t—s 
a 





__texp(—2tkf) (eset eminh eta) 2f coe Hn) 4, 


7 (t—s)?+4f? 


—b 


=m,(t) (—a <t <a) 


gz(t)— & (eden k(t—s) ag 


and 





~b 
__texp(—2ikf) [ g,(s)[(t—s)sinh k(t—s) + 2if cosh k(t—s)| ra 
- j (a af? 





=m,{t) (—b <t <b). 


The functions m,(t) and m,(t) are known continuous even functions given 


by 
dw (n= 1,2 





m,(t) = = 25 | af,,(a)cosh k(t?—w)! 


(?@—w*)t 
where m,(t) is defined for 0 < ¢ <a and m,(t) for 0 <t < b, and 
m,(t) = m,(—t)(—a <t<0), —_m,(t) = m,(—t) (—b <t < 0). 


The method given for two perfectly rigid or soft disks can readily be 
extended to diffraction problems for three or more disks. The Helmholtz 
function V(@,z) for such a problem is taken as a sum of contour integrals, 
each integral being associated with one of the disks and representing in the 
forms (3.5) or (3.4) the Helmholtz function for a fictitious distribution of 
doublets or sources on the disk. The set of Fredholm equations governing 


the problem can then be obtained by means of the boundary conditions 
on the disks. 
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SUMMARY 

Approximate expressions for the natural modes and frequencies of a uniform 
inextensible chain suspended from two points at the same level were obtained in 
a recent paper by Saxon and Cahn (1), who based their work on the solution of the 
fourth-order differential equation for the tangential displacement v. In the present 
paper the writer has obtained equivalent expressions using the angular displace- 
ment @ as independent variable. The differential equation for @ is of the second 
order, and was originally obtained by Routh (3) in a more general form for a non- 
uniform chain at least seventy years ago. 


1. Notation 

Ox, Oy are rectangular axes, Ox being the directrix of the catenary 
assumed by the chain when at rest. : 

L = total length of chain. 
m == mass of chain per unit length. 
T = tension at any point P in the equilibrium position. 
H = horizontal component of 7’. 

8 = are of chain measured from its mid-point, positive towards B. 

ys = angle between tangent at P and the axis Ox when the chain is at 

rest, with ¢ = +a at B, A.t 


The following relations hold for the catenary: 


L? 


8s = ctany#, where c = parameter of catenary = 4L cota = 


b = central dip of chain = c(sec a—1). 
y = csecy = ccosha/c. 
T = Hsecy, H = mge. 


Let the chain be in motion in a principal mode of frequency N (= w/2z7), 
and let 2+£sinw?t, y+nsinwt be the coordinates of P at time ft, 
%+@sinwi = angle made by the tangent at P with the axis Oz, 
T+ U sin wt = tension at P, and let U = GH sec yw. 

The motion of a string or chain was discussed very thoroughly by 
Routh in his classic Advanced Rigid Dynamics (3). In the chapter on this 

+ Saxon and Cahn used the symbol a for this quantity, and a, for its maximum value. 


[Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 1, 1961] 
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subject he defines the motion in terms of horizontal and vertical displace- 
ments, and this procedure is followed by Saxon and Cahn, who then 
derive the corresponding equations for the normal and tangential dis- 
placements (u,v). Routh, however, had the habit of putting many of his 
most interesting applications of general theory in small print, and in 
section 604, half-hidden in this manner, he develops the equation of motion 
of a ‘heavy slack heterogeneous chain’ in its intrinsic form.t The govern- 
ing differential equation is only of the second order in this case, and in 
Routh’s notation takes the form 


hp geosa/ée*d 
= — +2C}, 10’)? 
where a = ¢, ¢ = Osinwt and p is the curvature of the chain at rest 
(= ds/d{) in the present writer’s notation. For motion in a principal 
mode, C = Asinwt, where A is an arbitrary constant. 


ty 





7 sin we 








Fic. 1 

Routh did not consider the general solution of this equation, even 
approximately, but went on to investigate the motion of a chain whose 
mass distribution is such that it hangs in equilibrium in the form of a 
cycloid, a problem which can be solved exactly. 

The present paper is concerned with the approximate solution of Routh’s 
equation when the chain is uniform. We then have p = csec*%, and 
equation (10’) immediately transforms to 
0 (2 ge0%s-+-4)0-+2A ==: @ (1) 
where A? = w*e/g. 


+ The fifth edition of Routh’s work was published in 1892, and the first edition thirty- 
two years earlier. Routh’s intrinsic equation of motion of a chain is therefore at least 
seventy years old. ¢ Equations in (3), section 604, are here denoted by primes. 
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2. Boundary conditions 


The ‘determination of the whole motion’ was discussed by Routh in 
section 605 of (3), but only in rather general terms, and for the present 
purpose it is necessary to give a brief outline of the derivation of the 
equations of motion, in order to show the connexion between the variables 
6, €, and ». 

If we consider the element PP’ of the chain as a short inextensible link, 
the movement of P’ relative to P due to a small rotation @ of the link 

ives 
. =! = —@siny, and = bcos. (2) 


The equations of motion of PP’ in the z- and y-directions are respectively 


ua 


moa (€ sin wt) = x 6 P+ U sin wt)cos(b+-6 sin wt)}, 


my sin wt) = S(T+U sin wt)sin(y+ 6 sin wt)}—mg. 


Neglecting squares and products of small quantities, these reduce to 
_ sat - £ - (G—Btany) 


_ mewn 
H 


Putting ds = csec*y dy and eliminating £ and y between equations (2) and 
(3), it is found that d*G/df? = 2(d0/ds), whence 
b 
G= 2 | 6 d+ A+B. (4) 
0 
This is essentially the same as Routh’s equation (9’). In the process of 
elimination we obtain two relations between G and @, and substituting 
equation (4) in either of these immediately leads to equation (1). The 
general solution of the latter equation will contain two arbitrary constants, 
in addition to A, and the expressions for and 7 will thus involve four 
arbitrary quantities. The vanishing of ¢ and » at each end of the chain 
then gives four homogeneous relations between these quantities; but since 
the chain is symmetrical about its mid-point the natural modes may be 
split up into symmetrical and anti-symmetrical types, each involving only 
two constants. The frequency equation then follows in the usual manner. 


(3) 
and 


> £(G tan %+ 6) 


3. The approximate solution of equation (1) 
Since the parameter A increases with the mode number a solution in 
descending powers of A is likely to be the most useful form. A particular 
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integral of the equation is evidently 
A cos*ys 
oe O(;3) 


while the complementary function may be obtained by discussing the 
general equation ae 


dif? 
Making the Riccati transformation @ = K exp| f d(w) ay}, equation (5) 
0 


+(A°F*(h)+-a?)0 = 0. (5) 


becomes 


$+ P2+MFHY) +a? = 0, (6) 
where primes denote derivatives with respect to y. For A large, 
@ = +iAF. 


Hence, we may write 


, ik 8 l 
$= LIOP+O4 +5 +0(;5} 


Substituting in equation (6), expanding in descending powers of A, and 
equating coefficients to zero, we find that 


Q@=-= R= Bm abe om sar) etc. 


oF” oF’ © = —aglarp ") 


Putting ¢ = B+ty, where B and y are real functions of ¥%, we find on 
substituting in equation (6) that B = —y’/2y, whence 


exrf | Bad] 7 


Further, put jra— = A a nig) +” +0(53) 


then 0 = Ky-te+tA (8) 
and 6 = Ky-(B-+iy)e**. (9) 


The real and imaginary parts of equation (8) provide two linearly inde- 
pendent real solutions of the differential equation. 


In the present case, F(ys) = sec! and a? = 4. Thus, 
Q@ = —}jtany, R = }cos'y(13—}tan%J), ete., 


B = —Htany-+0( 55} 
on (13 — ?} tan*) +015 i) 


y = Asecly + 
v7 


p(y) = [ sectpdy, and gf) =} j (13—} tan%)ooslys di. 


ft) 
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2. Boundary conditions 

The ‘determination of the whole motion’ was discussed by Routh in 
section 605 of (3), but only in rather general terms, and for the present 
purpose it is necessary to give a brief outline of the derivation of the 
equations of motion, in order to show the connexion between the variables 
6, €, and ». 

If we consider the element PP’ of the chain as a short inextensible link, 
the movement of P’ relative to P due to a small rotation @ of the link 

ves 
” 3 = —@siny, and . = Ocosy. 


The equations of motion of PP’ in the z- and y-directions are respectively 


 . 7 oy . 
m—, (€sin wt) = ag (Fr +l sin wt)cos(+-@ sin wt)}, 


a2 
m= (7 sin wt) = £(T+U sin wt)sin(b+6 sin wt)}—mg. 


Neglecting squares and products of small quantities, these reduce to 
on = 


= gga tang) 


2 
and a = 4 (Gtany+6) 
Putting ds = csec*y dy and eliminating £ and » between equations (2) and 
(3), it is found that d?G/dy? = 2(d6/ds), whence 
ob 
G=2| Odp+Ad+B. (4) 
0 


(3) 


This is essentially the same as Routh’s equation (9’). In the process of 
elimination we obtain two relations between G and @, and substituting 
equation (4) in either of these immediately leads to equation (1). The 
general solution of the latter equation will contain two arbitrary constants, 
in addition to A, and the expressions for € and 7 will thus involve four 
arbitrary quantities. The vanishing of ¢ and 7 at each end of the chain 
then gives four homogeneous relations between these quantities; but since 
the chain is symmetrical about its mid-point the natural modes may be 
split up into symmetrical and anti-symmetrical types, each involving only 
two constants. The frequency equation then follows in the usual manner. 


3. The approximate solution of equation (1) 
Since the parameter A increases with the mode number a solution in 
descending powers of A is likely to be the most useful form. A particular 
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integral of the equation is evidently 
es ae co 4 Ol; 4} 


while the complementary function may be obtained by discussing the 
general equation Pry) 


dy? 
Making the Riccati transformation @ = K exp| j () ay, equation (5) 


+ (A2F?(y)+a*)? = 0. (5) 


becomes 


¢'+¢?+A2F2(4)+a? = 0, (6) 
where primes denote derivatives with respect to y. For A large, 


$m +iAF. 


Hence, we may write 


$= +iF+Qs'" ++) 


Substituting in equation (6), ee in descending powers of A, and 
equating coefficients to zero, we find that 


6, 2 8 Se Bs slop) ete. 


oF’ oF’ = —ay\ar)’ ( 


Putting ¢ = B+ty, where f and y are real functions of %, we find on 
substituting in equation (6) that B = —+y’/2y, whence 


ef fas] cy 

ub 
Further, put | ydf=A= rw) +42) + 0(;5), 
then 6 = Ky-te+#4 (8) 
and 0 = Ky-*(Btiy)e*. (9) 
The real and imaginary parts of equation (8) provide two linearly inde- 


pendent real solutions of the differential equation. 
In the present case, F(y) = sec!ys and a? = 4. Thus, 


Q = —jtany, R = }cos!{(13—} tan), etce., 
p= —ttany+0(5) 


y = Asecly + 2°54 (13 ? tan*y)+ O55): 
ob 
Pip) = [ seclbdy, and gi) =} (13—} tan%p)oosly dy. 


0 0 





122 W. J. GOODEY 
The integrals p and q may be expressed in terms of elliptic integrals and 


¢ 
trigonometrical functions (see Appendix). The integral f 6d, of the 


0 
complementary function may be found by writing it in the form 6T'+ a 
constant. Then ’+¢r =1. (10) 


For A large, ¢ = +iAF, so that T = +Fi1/AF. Hence, we may write 
r) r 1 ’ 


We then find from equation (10) that r= —3F’/2F°, so that in the 
present case r = —{sinycos*%. Equation (4) now becomes 
G = 2Kv-*(lLim)e*#4 4 Ap+iB. (12) 
The arbitrary constant is here written in the form +iB to make the real 
part of G an odd function of % and the imaginary part an even function. 
(This only applies where, as in the present case, A is an odd function of #.) 
In the odd modes, @ is an even function of %, and G is quite clearly an 
odd function. The odd modes are then given by the real parts of equations 
(8), (9), and (12), and the even modes by the imaginary parts. 


4. Odd modes 
In the fundamental and higher odd modes the mid-point of the chain 

moves almost horizontally in the vertical plane through the end-points, 
like the bob of a simpie pendulum; @ is therefore an even function of ¥, 
and the motion is defined by the real parts of equations (8), (9), and (12), 
to which must be added, in the case of (8) and (9), terms due to the 
particular integral. Following Saxon and Cahn, we reject second and 
higher powers of 1/A, whence 

6 = Ky-*tcosA 

0’ = Ky-*(B cos A—ysin A) 
and G = —2Ky-msin A+ Ad 
while, from (4), G’ = 2Ky-tcosA+A 


(13) 


Since € and » vanish at each end of the chain, we must have, from 
equations (3), 

7’ —6' tan y—Osec*y = G’ tand+Gsecf+6’ = 0, when ¢ = +a. 
Using equations (13), the condition that these last two equations are 
consistent gives 

(J—¢ tan*)— tan 4—B—By tan %)cos A+ 
+(2m+y+ytany¢)sinA=0 (f= +a). 
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The term 2m may now be omitted, since m = O(y/A*). Finally, putting 
B = —}tana, y = Asecla, and A = Ap,+q,/A, this result may be put in 
the form (cf. equation (29) of (1)) 


A. ie SE ') | 4 
tan(ap,+%3) = iseecpeoes } sin acosta} +0 Wh (14) 


Writing (« cos'a)(a sin a+cos «)-!—} sin a costa = p,, 


and observing that when A is large Ap, must be near nz, equation (14) 
may be solved approximately in descending powers of n. We thence find 


Al Nan a(y — dart PatPs + 0(-,)| (n = 1,2,...). (15) 
Px n® 3 n4}| 
This may be compared with equation (31) of (1). The quantity p, is the 
same as Saxon and Cahn’s f, but in place of (¢,+-p,) they have (g+-h). 
However, _ 
(da t+Pa)—(g+h) = } ( (2—tan*s)cos! d—} sin a costa = 0, 
0 
so that our results are in agreement. 
When the central dip of the chain is small compared with its length, 
a simple approximation for the frequencies of the odd modes may be 
found by expanding the right-hand side of equation (15) in ascending 
powers of a, ignoring fourth and higher powers, and finally putting 
x = 4b/L. Since N = (A/2z)(g/c)*, and c = b(seca—1)-', we have 


>, ® (g\h{ 38  2\b?) P 
W we (2) — (35+ 3) za (18) 
When n = 1, this becomes 


Ve (5) (1-482 73) (17) 


2v2\b | Fas 


This result may be compared with Pugsley’s empirical estimate of 


> 1 /g\t b? 
we, nals) (1-375) 


in (2), equation (21), and also with equation (15) of the same paper, giving 


ee (F)*(: oe B) 
2V2\b 7D 


for a cycloidal chain of small central dip. Pugsley’s formula does in fact 
show better agreement with the frequency calculated from equation (15) 
when 6/L is not very small, and appears to be a very good approximation 
up to about a = 60° (see Fig. 2). 
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5. Even modes 

In the even modes the mid-point of the chain moves up and down, 
6 being an odd function of %. The motion is therefore defined by the 
imaginary parts of equations (8), (9), and (12), and there is no particular 
integral, since A = 0. 

The frequency equation is now very simply obtained from the first of 


0-6 ‘ 
wey" 


0-5} 
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Fic. 2 

equations (3) by putting = 0 when ¢ = +a. The condition that » = 0 
when % = +a serves only to determine the ratio B/K. Thus we now 
mare @ = Ky-tsinA = 36’, 
and # = Ky~*(Bsin A+ cos A), 
whence 

Ky~*(2sin A—f tan sin A—+y tan ¢ cos A—sec*ys sin A) = 0, 
when ¢ = +a, leading to 


cot( Ap, + 3) = Ceres + 0(;5) (18) 





A Asin « 8 


Writing (1—} tan*a)cos'a cosec a = o,, and observing that when A is large 
Ap,, must be near (n-+-})7, equation (18) yields the approximate solution 
(n+4)r{,  (atoa)p l 
id Jy _ 4a wha Of. , 
; Pa | (ntan® t (nif 
There is no solution corresponding to n = 0, for since the chain is in- 


extensible the quantity @, and therefore sin A, must pass through zero at 
least once in the range 0 < yf < a. 





= 1,2....). (19) 
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Equation (19) above corresponds to equation (32) of (1). Saxon and 
Cahn’s expression (g+&) is slightly different in form from (q¢,+-¢,), but 
these quantities are in fact equal, for the same reason as in the case of 
the odd modes. 

When the central dip of the chain is small compared with its length 
a simple approximation for the frequencies of the even modes may be 
found from equation (19) by expanding in ascending powers of «, as was 
done for the odd modes. Thus it is found that 


When n = 1, this becomes 


N = a! } 
2v2 \b 


This may be compared with Pugsley’s estimate of 


,_, L4/(g\h 1-5b? 
Ne 3al5) (!— ae) 


in (2), equation (22), and also with equation (16) of the same paper, giving 
yw 1°43/g\4/, _ 1-58662 
~ ™ a2\b PB 
for a cycloidal chain of small central dip. In this case equation (21) appears 
to be a satisfactory approximation to the frequency of the second mode 


up to about a = 60°, being slightly better than Pugsley’s approximation 
(see Fig. 2). 


6. The oscillations of a vertical chain 

If « = 4n, the chain hangs in a vertical loop, and the odd modes will 
clearly be the same as those of a single piece of chain of length 4 
suspended from one end, since the two halves of the loop swing together. 
The even modes will be the same as the modes of a single piece of chain 
suspended from one end, with the other end constrained to remain 
vertically below the point of suspension. For a principal mode of frequency 
w/2m, the sideways motion of any point of the chain is £sin wt, and since 
the chain is inextensible there is negligible vertical motion, and the tension 
in the chain therefore differs inappreciably from the static tension mgs. 
The equation of sideways motion then takes the form 


m= (€sin wt) = Ms mgs © (€ sin wt)|, 
ot? és\ os j 


oda) to = @ 
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Putting s = gz*/4w*, this equation is transformed to 
dé i1dé 
antiga tt= © 

which is Bessel’s equation of order zero. 

This investigation is essentially the same as that given by Lamb (5) for 
a single piece of chain, and our purpose here is to relate it to the behaviour 
of a vertical loop. 

In the notation of (4), the general solution of equation (22) is 

§ = Ad(2wy/(8/g)) + BNg(204(6/9))- 
But N, > as 8-0, whence B= 0. The condition that ¢ = 0 when 
8 = 4L = 6 thus becomes J,{2w,/(b/g)} = 0, giving 
N= rad x zeros of Jp. (24) 

The frequency of the fundamental mode, in which the two halves of the 
loop swing together, is approximately 0-191(g/b)t. This is also, however, 
the frequency of the first even mode (i.e. the second mode of the loop), 
since by taking B infinitesimally small we can make £ = 0 at s = 0. 
This will have no effect on the frequency, since ¢ is still equal to 
AJ,j\2w,/(8/g)} except at s = 0. The distinction between even and odd 
modes has now disappeared, and the physical explanation appears to be 
that, since the tension in the chain is zero at the bottom of the loop, both 
halves of the chain can swing independently of each other as long as the 
amplitude remains small, continuity of the loop being preserved in the 
even modes by the formation of a small right-angled bend at the lower 
end. The foregoing investigation provides a limiting point for the curves 
given in Fig. 2, and serves to indicate how far the approximations obtained 
from equations (15) and (19) may be taken. It appears that they are 
reasonably accurate at least up to a = 70°, as was suggested by Saxon and 
Cahn. 


APPENDIX 


‘ “ 
Note on the integrals p(y) = { sec) dy and g(x) = 4 f (13—} tan%))cos'y dy 
0 0 


These integrals may be expressed in terms of trigonometrical functions and the 
standard elliptic integrals of the first and second kinds. In the notation of (4), the 
latter are respectively 

a“ “ 
F(k,p) = | (1—k*sinty)-t dy and E(k,p) = | (1—k*sin*y)t du. 
6 6 
Putting cosy = cos*u, we have 


dp(p) = 2 8ec*u(1-+cos*u)-t dy. 
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{tan y./(1+cos*u)} = (sec*z-+cos*z)(1+cos*x)-+ 
ld : l E 
=5 3 +~2(1—}sin*u)!— yp (l—dsinty)-t. 


a 
du 


1/1 1 
ence = i 2 eee Si. at — oe 3 
Hence, p = 2tanpy(1-+cos 1+ P(e) v2 B{ =H) (Al) 


v 


(In the tables on pp. 134-9 and 141-4 of (4), the parameter k is replaced by 
a are sin k; thus in the present case a = 45°.) 


v © 
To find g, put f cos’ dp = J, and f secty dy = J,, whence 
a d 


v 
I,—I, = | tan*pcos'y) dip. 
6 
Also dl, = 2cos*y(1+cos*u)-t du and di, = 2(1+cos*u)-tdu. Hence 


{1 
I, = v2 F( +n). 
Further, 
3dh,_1 di, 


| FY : 
r= {sin 2 cos j./(1-+cos*4)} = (3 cos*~— 1)(1+ cos*u)~+ 3 du ~~ du’ 


, 2 ./ 1 
giving I, = $sin pcos p/(1+-cos*z) 4 5 F(—.n)- 


Finally, q = #4(551,—3/,). 
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CORRIGENDUM 


MACK: ‘ROUTH TEST FUNCTION METHODS FOR THE 
NUMERICAL SOLUTION OF POLYNOMIAL EQUATIONS 


[Received 30 April 1960] 


My attention has been drawn to the paragraph starting ‘Finally, it is 
worth ...’ on page 370 of my paper (1) ‘Routh Test Function Methods for 
the Numerical Solution of Polynomial Equations’, and I regret that this 
conveys the impression that convergence cannot be guaranteed with any 
of the test functions given by Frazer and Duncan (2). An amended para- 
graph now follows: 

Finally it is worth pointing out the advantage of the x-method over 
Frazer and Duncan’s test-function methods (2). They give three test 
functions, namely, 7',_,, which is a determinant of order n—1, and two 
alternative functions, which are easier to compute, 7",_, and 7%_,. In 
the notation of the present paper, 

T= CP CPCPH ... Ce, 
a form which is also very convenient for computation. Now T7,,_, is 
polynomial in z and if it has different signs for x = x, and x = 2, oe 
there must be an odd number of factors of f(z) of the form z?—az-+-b where 
the a’s lie between 2z, and 22x,, the corresponding zeros being complex or 
not according as a? $ 4b. In such cases inverse interpolation for x from 
T,,., will produce rapid convergence in the determination of a correct 
factor z2—az+6. However, the function T’,_, suffers from the same 
disadvantage as the rational function C\*-;” in that both can change sign 
through an infinity. The function 7',_, “Ben an extraneous polynomial 
factor in x which can also change sign. Hence, in general, the use of one 
of the functions 7%,_,, T,_, or C&5" will not guarantee convergence until 
a close first approximation has “a obtained unless a more detailed in- 
vestigation establishing the true nature of any detected change of sign 
is also carried out. All the test functions have one common disadvantage, 
namely, that if there is an even number of factors z?—az+-b, where the a’s 
lie in (2a,, 2a,), the sign of the functions may not differ for x = x, and 
x = 2,. With the x-method the existence or non-existence of zeros with 
real parts in (x,,x,) is clearly demonstrated by computing the two cases 
x = 2, and z = 23. 
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not according as a? = 4b. In such cases inverse interpolation for x from 
T|,_, will produce rapid convergence in the determination of a correct 
factor z?—az+6. However, the function 7%,_, suffers from the same 
disadvantage as the rational function C’*—" in that both can change sign 
through an infinity. The function 7",_, has an extraneous polynomial 
factor in x which can also change sign. Hence, in general, the use of one 
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namely, that if there is an even number of factors z?—az-+-b, where the a’s 
lie in (2%,, 2a,), the sign of the functions may not differ for x = x, and 
x = 2. With the x-method the existence or non-existence of zeros with 
real parts in (x,,2,) is clearly demonstrated by computing the two cases 
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